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I.  Introduction 


With  the  development  of  near-infrared  (NIR)  laser  diodes,  the 
synthesis  of  fluorescent  dyes  with  excitation  and  emission  spectra  in  the 
NIR  wavelength  regime  has  accelerated  in  the  past  decade.  Examples  of 
fluorescent  dye  development  include  (i)  microscopy  dyes  for  sensing 
cellular  metabolic  and  structural  constituents;  (ii)  DNA-binding  dyes  for 
automated  DNA  sequencing;  and  (iii)  photodynamic  agents  for  anti¬ 
cancer,  therapeutic  action.  Armed  with  the  knowledge  that  NIR  light  can 
penetrate  several  centimeters  into  tissue  [for  review  see  reference  (1)], 
one  can  speculate  that  these  and  other  designer  fluorescent  dyes  may 
offer  diagnostic  opportunities  when  used  in  conjunction  with  optical 
imaging  techniques  as  contrast  agents  for  the  detection  of  diseased 
tissues,  specifically  breast  cancer. 

The  detection  of  disease  via  the  use  of  fluorescent  dyes  and 
photodynamic  drugs  as  contrast  agents  for  the  detection  of  disease  has 
been  proposed  by  several  investigators.  However,  a  long-standing 
problem  has  been  the  low  uptake  or  "leakage"  into  diseased  tissues 
providing  insufficient  contrast  for  detection  of  disease  -  making  it 
difficult  to  detect  the  diseased  tissue  signal  from  the  background  tissue 
signal.  While  targeted  delivery  may  improve  uptake  ratios  and  contrast 
imparted  by  these  agents,  drug  and  contrast  agent  targeting  approaches 
have  been  rather  elusive  even  in  conventional  imaging  modalities  such  as 
magnetic  resonance  imaging  and  x-ray  computed  tomography. 
Investigators  have  sought  to  improve  optical  contrast  independent  of 
uptake  ratio  by  employing  dyes  which  fluoresce  differently  in  diseased 
and  normal  tissues.  The  use  of  agents  whose  emission  characteristics 
vary  with  tissue  pH  (2,  3)  and  02  (4)  may  not  only  enhance  detection  of 
diseased  tissues  by  the  nature  of  differing  fluorescent  properties,  but 
may  also  contain  diagnostic  information  regarding  the  disease. 

The  difficulty,  however,  lies  in  using  the  multiply  scattered  light 
detected  at  the  tissue-air  interface  to  reconstruct  an  image  which 
differentiates  internal  tissues  on  the  basis  of  uptake  and  the  fluorescent 
properties  of  quantum  yield,  and  lifetime.  In  ongoing  work  in  this  and 
other  laboratories,  investigators  have  demonstrated  the  capacity  for 
fluorescent  lifetime  and  quantum  yield  imaging  from  frequency-domain 
measurements  of  modulation  phase-shift  and  amplitude  demodulation 
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conducted  at  the  surface  of  a  simulated  tissue  phantom  (5-7).  Under 
this  support,  the  aims  of  this  two  year  exploratory  research  project  are: 

(1)  To  develop  inverse  imaging  algorithms  providing  lifetime  and  yield 
mapping  of  tissues  as  thick  as  five  centimeters; 

(2)  To  measure  the  temporal,  re-emitted  excitation  and  fluorescence  from 
tissue  phantoms  using  frequency-domain  techniques;  and 

(3)  To  use  these  measurements  as  input  into  the  inverse  imaging 
algorithm  for  reconstruction  of  lifetime  and  yield  maps. 

Finally,  in  this  current  last  year  of  the  project,  we  will  accomplish  the 
last  and  final  aim  of  the  project, 

(4)  To  test  computational  and  experimental  approaches  using  dyes  which 
enable  analyte  concentration  determination  from  measurement  of 
lifetime  and  yield  within  a  tissue-mimicking  phantom. 

Specifically,  we  seek  to  develop  fluorescence  lifetime  imaging 
techniques  for  screening  in  a  population  of  women  who  are  at  high  risk 
for  the  disease,  and  for  those  women  who  require  prognostic  indicators  of 
lymphatic  involvement.  Under  other  support,  the  development  of 
imaging  instrumentation  and  contrast  agents  for  identification  of  disease 
is  ongoing.  This  DOD  support  focuses  upon  the  development  of  inverse 
solution  procedures  and  coupling  of  measurements  to  image  algorithms 
for  the  detection  of  disease.  These  specific  aims  are  the  foundations  of 
the  technology. 

In  the  following,  we  outline  the  computational  and  experimental 
methods  as  well  as  the  results  obtained  in  this  first  year  of  the 
exploratory  research  grant. 

II.  Computational  methods 

During  the  past  year,  computational  efforts  for  reconstructing 
lifetime  and  quantum  yield  from  measurements  of  phase  and  amplitude 
modulation  have  been  spearheaded  by  Dilip  Y.  Paithankar,  and  more 
recently  Tamara  L.  Troy  in  the  Photon  Migration  Laboratory.  Both 
studies,  although  separate,  iteratively  solved  the  diffusion  equations 
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describing  light  propagation  and  fluorescence  generation  in  a  random 
medium  as  the  spatial  distribution  of  lifetime  and  quantum  efficiency 
were  altered.  The  map  of  these  fluorescent  optical  properties  that  best 
matched  simulated  measurements  and  predictions  from  the  diffusion 
equations  provided  the  optical  image  of  interior  heterogeneities  from 
exterior  measurements.  The  preliminary  work  described  in  the 
publication,  “Imaging  of  fluorescent  yield  and  lifetime  from  multiply 
scattered  light  re-emitted  from  random  media,”  involved  reconstruction  of 
lifetime,  x,  and  the  product  of  the  quantum  efficiency  and  absorption 
coefficient  owing  to  fluorophore,  >  within  a  circular  2-D  phantom 

with  20  detectors  and  4  sources  of  modulated  light.  This  latter  quantity, 
(l)(iax  m  >  is  related  to  the  local  concentration  of  fluorophore,  [C],  via  the 

extinction  coefficient,  s: 


=MC]  (i) 

and  convolutes  changes  in  contrast  agent  uptake  with  changes  in 
quantum  yield  due  to  differing  biochemical  environments.  Details  of  the 
computational  approach  and  results  are  shown  in  the  publication 
included  in  the  Appendix.  In  brief,  the  results  show  that  it  is  possible  to 
reconstruct  interior  images  of  lifetime,  x,  and  fluorescent  yield,  4>Hax_>m  > 

from  exterior  measurements  of  phase-shift  and  amplitude  modulation 
made  at  the  tissue-air  interface. 

However,  uptake  or  local  concentration  of  fluorophore,  [C],  may  not 
be  specific  to  disease.  The  biochemistry  probed  quantum  yield,  <j>, 
changes  of  designer  fluorophores  may  very  well  be.  Consequently,  to 
enhance  the  specificity  for  detection  of  breast  cancer  in  women  who  are 
at  high  risk,  or  who  show  positive  mammograms,  we  sought  to  develop 
an  imaging  technique  that  focused  upon  specificity  to  disease  as  detected 
by  <|)  rather  than  the  product,  4>hax_^m  .  More  recently,  Troy  and  Sevick- 

Muraca,  developed  an  inversion  technique  that  could  segregate  changes 
in  lifetime  and  quantum  efficiency  from  differences  in  local 
concentrations  of  fluorescent  contrast  agents.  This  is  especially 
important  for  the  specific  identification  of  disease  based  upon  local 
biochemical  environment. 
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In  this  second  approach  in  the  laboratory,  the  uptake  of 
fluorophore  is  independently  assessed  by  reconstruction  of  absorption 
owing  to  fluorophore  alone,  or  pa  m  •  The  reconstruction  is  performed 

by  simulating  measurements  of  phase-shift  and  amplitude  modulation  at 
excitation  wavelength  prior  to  and  following  administration  of  a  contrast 
agent.  With  the  exterior  values  of  Pax_>m  determined,  simulated 

measurements  of  phase-shift  and  amplitude  modulation  at  the  emission 
wavelengths  are  then  used  to  reconstruct  values  of  lifetime,  x,  and 
quantum  efficiency,  <j>. 

In  the  most  recent  work  by  Sevick  and  Troy,  the  simulated 
phantom  mimics  the  experimental  system  in  the  Photon  Migration 
Laboratory  consisting  of  a  square  rectangular  phantom  whose  surface  is 
imaged  onto  the  photocathode  of  a  gain-modulated  image  intensified 
CCD  camera.  While  measurements  of  fluorescent  and  excitation  phase- 
shift  and  amplitude  modulation  across  one  face  of  the  phantom  are 
possible  at  a  resolution  of  512  x  512  (8,9);  we  simulate  the  use  of  1x17 
detection  pixels  in  our  inversions.  With  four  sides  imaged,  there  are  56 
(4  x  17  -  2)  detectors  positioned  around  the  periphery  and  4  sources 
located  midplane  on  each  of  the  four  sides. 

Under  conditions  in  which  contrast  is  obtained  from  a  10:1  uptake 
alone,  i.e.,  pa  >  the  2-D  detection  of  a  0.06  cm2  heterogeneity 

embedded  in  a  16  cm2  phantom  can  be  distinguished  from  a  somewhat 
noisy  background,  as  shown  in  Figure  1.  This  uptake  ratio  is  the  best 
physiologic  uptake  reported  for  a  photodynamic  agent.  It  is  important  to 
note  that  the  data  used  in  the  reconstruction  of  Figure  1  is  simulated 
and  not  experimental.  The  data  are  simulated  with  an  ample  noise  level 
of  0.1  degree  in  phase-shift  and  1%  in  amplitude  demodulation  at  the 
excitation  wavelength  at  100  MHz  modulation  frequency. 

On  the  otherhand,  when  additional  contrast  is  provided  by  a 
decrease  in  lifetime  from  10  nanoseconds  in  the  surrounding  normal 
tissues,  to  1  nanosecond  in  the  simulated  diseased  heterogeneity,  the 
contrast  is  enhanced  providing  better  detection  of  the  heterogeneity 
(Figure  2).  This  latter  result  assumes  that  the  uptake  ratio  has  been 
correctly  detected  from  excitation  measurements.  It  is  important  to  note 
that  fluorescent  measurements  are  simulated  with  an  ample  noise  level 
of  0.1  degree  in  phase-shift  and  1%  in  amplitude  demodulation. 
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Less  encouraging,  but  still  promising  are  the  results  on  quantum 
efficiency.  Albeit  a  noisy  image  reconstruction,  Figure  3  illustrates  the 
additional  contrast  provided  by  an  increase  in  quantum  yield  from  a 
value  of  0.01  in  the  surrounding  normal  tissues,  to  a  value  of  0.1  in  the 
diseased  tissue.  These  results  are  again  provided  for  a  physiological 
uptake  ratio  of  10:1.  It  is  important  to  note  that  DNA-binding  dyes, 
although  not  candidates  for  contrast  agent  administration  in  patients, 
have  a  quantum  yield  change  of  one  thousand  fold  difference  in  healthy 
and  dying  cells. 

These  results  are  encouraging,  suggestive  that  fluorescence  lifetime 
and  quantum  efficiency  can  be  exploited  to  provide  contrast  beyond  that 
provided  by  uptake  changes. 

In  year  two,  work  continues  to  expand  imaging  algorithms  for  use 
of  multi-frequency,  multi-wavelength  measurements.  In  addition,  we  are 
exploring  the  use  of  non-iterative  methods  to  provide  rapid  image 
formation  for  real  time  clinical  use.  Our  work  in  this  domain  is  not 
described  here,  but  rather  in  the  manuscript  accepted  for  publication 
and  found  in  the  Appendix. 

III.  Experimental  Methods  and  Approaches 

In  order  to  experimentally  assess  the  contrast  provided  by 
fluorescence  and  absorption,  frequency-domain  measurements  were 
conducted  in  a  tissue-simulating  phantom  using  micromolar 
concentrations  of  IR-125  (or  indocyanine  green,  ICG)  (ACROS  Organics, 
Fairlawn,  NJ),  and  DTTCI  (ACROS  Organics,  Fairlawn,  NJ).  The  tissue 
phantom,  illustrated  in  Figure  4,  consisted  of  a  cylindrical  Plexiglas 
container,  20  cm  in  diameter,  30.5  cm  in  height,  and  filled  with  0.5% 
Intralipid  solution  (Kabi  Pharmacia,  Clayton,  NC)  to  mimic  tissue 
scattering  properties.  The  heterogeneity  consisted  of  a  9  mm  inner 
diameter  cylindrical  glass  container  in  which  the  same  scattering 
solution  was  held.  The  appropriate  concentrations  of  dye  were  added  to 
both  solutions  to  explore  the  issues  of  fluorescent  contrast  for  frequency- 
domain  photon  migration  imaging.  Measurements  of  phase-shift  and 
amplitude  demodulation  were  conducted  as  the  position  of  the  tissue 
mimicking  heterogeneity  was  precisely  moved  using  an  x-y  translation 
stage  (PMC200-P,  Newport,  Irvine,  CA)  to  various  distances  midpoint 
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between  the  source  and  detector.  Contrast  was  measured  as  a  function 
of  distance,  x,  between  the  inner  circumference  of  the  phantom  and  the 
leading  edge  of  the  heterogeneity  as  shown  in  Figure  4.  An  object 
position  of  zero  denotes  contact  between  the  phantom  wall  and  the 
heterogeneity.  Measurements  were  conducted  in  the  absence  of  dyes  to 
assess  the  contribution  of  the  vessel  on  measurements  of  phase-shift  and 
amplitude  modulation.  Measurements  of  excitation  light  showed  less 
than  a  one  degree  phase-shift  and  a  0.1%  modulation  change  owing  to 
the  presence  of  the  glass  vessel  filled  with  the  same  Intralipid  solution  as 
contained  in  the  surrounding  medium.  These  deviations  due  to  the 
vessel  wall  of  the  heterogeneity  occurred  only  when  it  was  close  to  of  the 
phantom  wall.  Since  these  differences  are  small  compared  to  the 
changes  owing  to  fluorescence  in  our  measurement,  we  did  not  account 
for  them  in  the  presentation  of  our  results. 

Indocyanine  green,  ICG  is  a  FDA  approved,  tricarbocyanine  dye 
used  typically  for  hepatic  function  studies  and  ophthalmic  angiography 
studies.  It  has  also  been  used  as  an  absorption  contrast  agent  for 
imaging  of  rat  gliomas  and  tumor  margins  (10)  and  as  a  fluorescent 
contrast  agent  for  the  estimation  of  burn  depth  (11).  Although  ICG  has 
absorption  and  fluorescence  peaks  at  764  nm  and  803  nm,  respectively, 
the  absorption  and  fluorescence  measurements  were  conducted  at  780 
nm  and  830  nm,  respectively,  to  be  consistent  with  the  characteristics  of 
DTTCI.  ICG  dosage  corresponds  to  uniform  systemic  concentrations  of  1- 
3  qM,  50-75  fold  lower  than  lethal  levels.  While  these  ICG 
concentrations  are  also  in  the  therapeutic  range  for  many  photodynamic 
agents,  we  have  found  that  the  quantum  efficiency  of  ICG  is  tenfold  or 
more  greater  than  many  photodynamic  agents.  Nonetheless,  in  the  work 
presented  herein,  we  assess  the  contrast  offered  by  ICG  for  photon 
migration  imaging  using  clinically  realistic  tissue  concentrations  of  ICG. 
Extrapolation  of  our  results  using  ICG  to  therapeutic  concentrations  of 
photodynamic  agents  must  account  for  the  differences  in  fluorescent 
yield. 

Measurements  of  phase-shift  and  amplitude  modulation  owing  to  dye 

Experimental  measurements  of  fluorescent  phase-shift  and 
amplitude  demodulation  were  conducted  at  modulation  frequencies  of 
80,  160,  and  240  MHz  using  an  ISS  (ISS,  Champaign,  IL)  based  system 
illustrated  in  Figure  5.  The  excitation  source  consisted  of  an  Argon- 
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pumped,  picosecond  pulsed  Ti:  Sapphire  laser  (Spectra- Physics, 
MountainView,  CA)  operating  at  780  nm,  at  a  laser  repetition  rate  of 
80.00000  MHz  and  an  average  power  of  1.3  W.  A  glass  microscope  slide 
was  used  to  split  the  beam,  directing  90%  of  the  light  into  a  1000  micron 
fiber  (HCP-M1000T-08,  Spectran  Specialty  Optics  Co.,  Avon,  CT)  coupled 
to  the  tissue  phantom  model.  The  remaining  10%  was  delivered  to  a 
reference  photomultiplier  tube  (PMT),  (R928,  Hamamatsu  Photonics)  to 
provide  a  reference  signal.  The  multiply  scattered  light  detected  at  the 
tissue  phantom  periphery  was  collected  using  a  second  1000  micron 
fiber  located  2.8  circumferential  centimeters  from  the  source  fiber.  The 
light  was  delivered  to  a  sampling  PMT  outfitted  with  either  neutral 
density  filters  (CVI  Laser  Corporation,  Albuquerque,  NM)  and  a  780+10 
nm  interference  filter  (model  #10-780-4-1.00,  CVI  Laser  Corporation, 
Albuquerque,  NM)  for  collecting  the  multiply  scattered,  excitation  light, 
or  an  830+10  nm  interference  filter  (model  #10-830-4- 1.00, CVI  Laser 
Corporation,  Albuquerque,  NM)  to  monitor  the  generated  fluorescence. 
The  PMTs  were  gain  modulated  at  a  harmonic  of  the  laser  repetition  rate 
plus  an  offset  frequency  of  100  Hz.  The  heterodyned  signal  was  then 
acquired  and  analyzed  for  phase-shift  and  amplitude  demodulation  using 
ISS  electronics  and  software. 

Measurement  of  extinction  coefficients ,  quantum  yields  and  lifetimes 

The  extinction  coefficients  of  DTTCI  and  ICG  were  determined  at 
780  nm  and  830  nm  using  a  double-beam  spectrophotometer  (Spectronic 
Genesys  5,  Milton  Roy,  Rochester,  NY)  from  a  series  of  micromolar 
samples  diluted  in  water.  Fluorescent  quantum  yield  measurements 
were  conducted  using  a  spectrofluorometer  (Fluorolog-2,  SPEX,  Edison, 
NJ).  These  values  were  determined  with  comparison  to  measurements  of 
ICG  in  DMSO  using  the  published  quantum  yield  at  780  and  830  nm 
excitation/ fluorescent  wavelengths  (12).  Lifetime  measurements  were 
conducted  with  phase  fluorimetry  using  IR-140  (ACROS  Organics, 
Fairlawn,  NJ)  as  a  standard.  Table  I  lists  the  extinction  coefficients, 
quantum  yield,  and  lifetime  of  the  two  dyes  in  water.  These  properties 
are  assumed  to  be  reflective  of  those  in  0.5%  Intralipid. 

Experimental  protocols 

Table  II  summarizes  the  experiments  involving  a  single  fluorescent 
heterogeneity  designed  to  evaluate  (i)  the  influence  of  lifetime  on 
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contrast,  (ii)  the  effect  of  uptake  ratio  on  the  fluorescent  and  absorption 
contrast,  and  (iii)  the  influence  of  fluorescent  surroundings  on  the 
contrast  enabling  detection  of  a  single  body. 

IV.  Experimental  Results  and  Discussion 

Influence  of  lifetime  on  fluorescent  contrast 

In  order  to  assess  the  influence  of  lifetime  on  contrast,  two 
Intralipid  solutions  containing  ICG  and  DTTCI  at  concentrations  of  2.0 
and  4.2  mM  respectively  were  employed  separately  as  heterogeneities 
within  the  tissue  phantom.  The  DTTCI  ad  ICG  concentrations  were 
chosen  since  they  reflected  similar  fluorescent  cross  sections  (<|>M.a  m  = 

4xl0-3  and  6xl0~3),  yet  different  lifetimes  (1.18  and  0.56  ns, 
respectively).  Figure  6a  and  6b  illustrate  the  phase-shift  difference  and 
amplitude  demodulation  for  the  excitation  and  emission  light  reported 
relative  to  the  absence  case  as  a  function  of  the  lateral  position  of  the  2.0 
mM  ICG  filled  heterogeneity.  The  closed  symbols  denote  fluorescent 
measurements  made  with  the  830  nm  interference  filter  while  the  open 
symbols  denote  excitation  measurements  conducted  with  neutral  density 
filters  and  the  780  nm  interference  filter.  Comparison  of  signals  shows 
that  greater  phase-shift  and  amplitude  demodulation  contrast  at 
distances  farther  from  the  periphery  is  afforded  by  fluorescence 
measurements.  This  occurs  due  to  the  lifetime  or  added  time-lag  of  the 
emission  process  that  is  inherent  in  fluorescence,  but  not  in  absorbance 
measurements.  In  addition,  the  fluorescent  contrast  owing  to  a 
heterogeneity  position  greater  than  4  cm  away  from  the  source  and 
detector  is  enhanced  at  lower  modulation  frequencies,  consistent  with 
the  observation  of  greater  excitation  photon  sampling  volume  at  lower 
modulation  frequencies  (12). 

Figures  7a  and  7b  contain  the  analogous  data  for  the  4.2  pM 
DTTCI  in  Intralipid  solution.  Since  similar  fluorescent  cross  sections  is 
available  in  the  two  measurements  involving  the  two  dyes,  the  difference 
in  contrast  between  the  DTTCI  and  ICG  fluorescence  measurements 
described  in  Figures  6  and  7  must  be  due  primarily  to  differences  in 
lifetime.  As  the  lifetime  of  the  fluorescent  dye  in  the  heterogeneity 
increases,  the  measured  amplitude  demodulation  should  decrease  and 
the  contrast  available  from  amplitude  demodulation  increases.  Indeed 
the  maximum  amplitude  demodulation  owing  to  DTTCI  is  approximately 
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0.5  while  for  ICG  is  0.7,  consistent  with  the  greater  lifetime  of  the  DTTCI 
dye.  In  addition,  as  the  lifetime  of  the  fluorescent  dye  in  the 
heterogeneity  increases  the  fluorescent  phase-shift  should  increase,  and 
the  contrast  available  from  phase-shift  increases.  Indeed  the  maximum 
contrast  afforded  by  phase-shift  owing  to  DTTCI  is  approximately  45 
degrees  while  for  ICG  is  approximately  30  degrees,  again  consistent  with 
the  greater  lifetime  of  DTTCI.  The  differences  in  absorption 
measurements  between  DTTCI  and  ICG  arise  due  to  the  differences  in 
absorption  cross  sections  of  the  two  solutions.  We  have  repeated  the 
experiments  shown  in  Figures  6  and  7  with  the  exception  that 
concentrations  yielding  identical  absorption  cross  sections  were 
employed.  The  results  show  similar  behavior  to  that  as  illustrated  in 
Figures  6  and  7  and  are  not  shown  here  for  brevity.  Nonetheless,  our 
results  show  that  frequency-domain  measurements  contain  information 
of  fluorescence  lifetime  and  can  be  used  to  provide  diagnostic 
information  if  lifetime  can  be  extracted  from  the  multiply  scattered  signal 
with  inverse  imaging  algorithms  as  proposed  by  Sevick-Muraca,  et  al.,  (8) 
and  O'Leary,  et  al.  (9). 

The  effect  of  uptake  on  fluorescent  and  absorption  contrast 

Since  the  administration  of  contrast  agents  can  hardly  expect  to 
result  in  perfect  uptake  into  tissue  regions  of  interest,  a  study  was 
conducted  to  investigate  the  reduction  of  contrast  owing  to  imperfect 
uptake.  Figures  8a  and  8b  illustrate  the  phase-shift  and  amplitude 
demodulation  as  a  function  of  position  of  a  heterogeneity  containing  a 
one  hundred  fold  greater  concentration  of  ICG  than  its  surroundings. 
Measurements  are  reported  as  a  function  of  the  heterogeneity  position 
away  from  the  source  and  detector.  The  concentration  of  ICG  in  the 
scattering  solution  comprising  the  heterogeneity  was  2.3  pM. 
Fluorescent  and  excitation  measurements  are  reported  relative  to  those 
conducted  in  the  absence  of  the  heterogeneity.  The  shaded  symbols 
denote  referenced  fluorescence  measurements  made  with  the  830  nm 
interference  filter  while  the  open  symbols  denote  referenced  excitation 
measurements  conducted  with  the  780  nm  interference  filter  and  neutral 
density  filters.  Upon  comparison  to  the  perfect  uptake  results  illustrated 
in  Figure  6,  one  finds  that  the  magnitude  of  contrast  is  dramatically 
reduced  when  the  partitioning  between  simulated  diseased  and  normal 
tissues  is  reduced  from  a  perfect  uptake  ratio  to  a  1:100  uptake  ratio 
from  the  surroundings  to  the  heterogeneity.  Nonetheless,  the  magnitude 
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of  contrast  offered  by  fluorescence  phase-shift  and  amplitude 
demodulation  measurements  still  exceeds  that  possible  by  absorption 
measurements  in  agreement  with  the  argument  presented  in  Section  2. 
This  again  can  be  attributed  to  the  additional  contrast  offered  by  the 
fluorescence  process.  In  addition,  the  experimental  measurements  of 
fluorescence  contrast  appears  to  be  more  sensitive  to  heterogeneity 
position  than  those  measurements  of  contrast  owing  only  to  absorption. 

Figures  9a  and  9b  illustrate  similar  measurements  when  the 
partitioning  between  the  heterogeneity  and  surrounding  scattering 
medium  is  tenfold.  Again,  the  absorption  measurements  indicate  the 
ability  to  detect  a  heterogeneity  located  between  zero  and  1  cm  from  the 
phantom  periphery  using  phase  and  amplitude  demodulation 
measurements,  while  fluorescent  measurements  afford  detection  of  the 
heterogeneity  located  twice  as  deep  within  the  phantom.  In  addition,  the 
sensitivity  of  phase  contrast  to  heterogeneity  position  appears  greater  for 
fluorescence  measurements  than  when  contrast  is  owing  to  absorption. 

Current  approaches  for  frequency-domain  optical  imaging  in 
tissues  involve  reconstructing  perturbations  in  tissue  absorption.  Since 
the  best  possible  contrast  that  can  be  offered  by  absorption  is  that  due  to 
a  perfect  absorber,  we  evaluated  the  contrast  due  to  fluorescence  in 
comparison  to  that  due  to  perfect  absorber.  Using  a  black  painted 
Plexiglas  rod  9  mm  in  diameter,  measurements  of  phase-shift  and 
amplitude  demodulation  of  the  multiply  scattered  excitation  light  at  780 
nm  were  conducted  as  the  perfect  absorber  was  moved  away  from  the 
source  and  detector.  Measurements  are  reported  relative  to  the 
"absence"  case.  Figure  10a  and  10b  illustrates  the  contrast  offered  by 
phase-shift  and  amplitude  demodulation  measurements  as  a  function  of 
heterogeneity  position  for  a  9  mm  diameter  perfect  absorber  and  a  9  mm 
diameter  tube  filled  with  a  3  pM  concentration  of  ICG.  In  the  case  of  the 
fluorescence  measurements,  the  phase-shift  and  amplitude 
demodulation  were  reported  relative  to  the  "absence"  case  in  which  the 
fluorescent  body  was  moved  more  than  5  cm  away  from  the  source  and 
detector.  These  results  point  to  the  superior  contrast  offered  by 
fluorescence  frequency-domain  approaches  since  the  contribution  of 
fluorescent  lifetime  is  included  in  the  measurement. 

Since  partitioning  of  contrast  agents  is  a  common  obstacle  to  even 
conventional  imaging  modalities,  one  might  expect  similar  difficulties 
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involving  fluorescent  contrast  agents  for  optical  imaging.  The  ability  to 
detect  a  fluorescent  body  in  scattering  medium  which  mimics  tissues  has 
been  demonstrated  by  several  investigators  However  no  studies  have 
addressed  the  contrast  available  when  there  is  imperfect  partitioning  in 
the  tissue  volume  of  interest.  We  have  demonstrated  that  the  contrast 
due  to  fluorescence  exceeds  that  due  to  absorption  even  under 
conditions  of  imperfect  partitioning.  The  advantage  for  fluorescent 
biomedical  optical  imaging  may  lie  in  the  additional  contrast  offered  by 
changing  the  quantum  yield  and  lifetime  of  contrast  agents  when  inside 
the  tissue  region  of  interest. 

Influence  of  lifetime  and  quantum  yield  in  heterogeneity  and 
surroundings 

In  order  to  evaluate  the  contribution  of  changing  lifetime  upon  the 
contrast  offered  by  fluorescent  frequency-domain  measurements, 
Intralipid  solutions  of  DTTCI  and  ICG  fluorescent  dyes  were  employed  in 
the  surroundings  and  the  heterogeneity  to  evaluate  the  impact  of 
changing  lifetimes  upon  measured  contrast.  In  all  experiments,  the 
concentrations  of  dyes  were  maintained  such  that  the  absorption  cross 
section  due  to  dye  in  the  heterogeneity  was  approximately  one  hundred 
times  that  in  the  surrounding  medium.  In  one  experiment,  the  DTTCI 
dye  solution  with  the  longer  lifetime  was  employed  in  the  heterogeneity 
while  the  ICG  solution  comprised  the  surrounding  medium.  In  the 
second  experiment,  the  long-lifetime  DTTCI  dye  solution  was  employed  as 
the  surrounding  medium  while  the  shorter  lifetime  ICG  solution  was 
used  in  the  heterogeneity. 

The  open  symbols  in  Figures  11  and  12  illustrate  the  phase-shift 
and  amplitude  demodulation  reported  at  80  MHz  relative  to  the  absence 
case  when  the  longer  lived  fluorescent  dye  was  partitioned  within  the 
heterogeneity  ([Cdttci]=4x10"6  M,  na  =1.7xl0'1  cm-1)  and  the  short¬ 
lived  dye  constituted  the  surroundings  ([Cicg]=1x10‘8  M,  Pax^m  =  1-3x10" 

3  cm-1).  The  filled  symbols  illustrate  similarly  obtained  data  with  the 
exception  that  the  longer  lived  fluorescent  dye  was  partitioned  within  the 
surroundings  ([Cdttci]=4x10‘8  M,  =l-7xl0'3  cm-1)  and  the  short¬ 

lived  dye  constituted  the  heterogeneity  ([Cicg]=1x10"6  M,  =1.3x10"1  cm* 
1).  The  filled  symbols  denote  referenced  fluorescent  measurements  made 
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with  the  830  nm  interference  filter  while  the  open  symbols  denote 
referenced  excitation  measurements  conducted  with  the  780  nm 
interference  filter.  There  are  several  points  to  note:  (1)  Since  the 
absorption  cross  sections  of  the  surroundings  and  heterogeneities 
remain  constant  between  the  two  experiments,  the  excitation  phase-shift 
and  amplitude  demodulation  measurement  in  Figures  1 1  and  12  show 
nearly  identical  phase  and  amplitude  demodulation  contrast;  (2)  The 
fluorescent  phase  contrast  exhibited  by  a  longer-lived,  fluorescent 
heterogeneity  in  a  shorter-lived,  fluorescent  surroundings  results  in  a 
greater  phase-shift  when  compared  to  the  fluorescent  phase  contrast 
exhibited  by  a  shorter  lived  fluorescent  heterogeneity  in  a  longer  lived 
fluorescent  surroundings  (Figure  11);  and  (3)  Even  though  the  quantum 
yield  of  DTTCI  is  three  times  and  the  lifetime  is  double  that  of  ICG  (see 
Table  I),  the  contrast  provided  by  amplitude  demodulation  results  in 
values  less  than  one  when  DTTCI  comprises  the  fluorescent  heterogeneity 
(Figure  12).  Conversely,  the  contrast  provided  when  ICG  comprises  the 
fluorescent  heterogeneity  results  in  amplitude  demodulation  values 
greater  than  one  (Figure  12). 

These  results  demonstrate  that  fluorescent  lifetime  and  quantum 
yield,  properties  that  are  dependent  upon  local  biochemical  environment, 
can  be  reflected  in  the  multiply  scattered  fluorescent  signals  detected 
from  the  periphery  of  tissue  simulating  phantoms.  Upon  application  of 
inverse  imaging  algorithms,  determinations  of  lifetime  and  yield  maps  of 
interior  heterogeneities  may  be  possible. 

V.  Summary 

The  performance  of  contrast  agents  for  conventional  imaging 
modalities  can  be  limited  by  the  less  than  perfect  partitioning  in  diseased 
tissues  of  interest.  The  use  of  fluorescent,  time-dependent  optical 
imaging  techniques  and  algorithms,  as  demonstrated  with  frequency- 
domain  approaches  shown  herein,  offers  the  ability  for  enhanced 
biomedical  contrast  owing  to  changing  fluorescent  properties  of  lifetime 
and  quantum  yield  when  localized  in  the  tissue  volume  of  interest.  In 
this  work,  we  have  pointed  out  the  possibility  for  enhanced  contrast 
provided  by  fluorescence  kinetics  and  suggest  that  contrast  agents 
designed  for  biomedical  optical  imaging  should  employ  these  concepts  for 
better  detection,  localization,  and  characterization  in  non-invasive  NIR 
imaging. 
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VI.  Year  Two  Activities 

In  year  2,  we  will  couple  the  inversion  algorithms  with 
experimental  measurements  conducted  with  a  image  intensifier  CCD 
camera.  Our  goal  will  be  to  perform  analyte  sensing  from  lifetime  and 
quantum  efficiency  measurements. 
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Dye  Xx  ^,m  e7g0  £830  $  x  (ns) 

(nm)  (nm)  (M*cm)~1  (M*cm)-1 _ 

780  830  43,000  5,500  0.034  1 

780  830  130,000  22,000  0.016  0 
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9het4ax.mhet 
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830 

2.0 

0 
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4x  1 0'3 

DTTCI  Perfect  Uptake 
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4.2 

0 

0.18 
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ICG  100:1  Uptake  Ratio 
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2.3 

0.023 

0.30 

5xl0-3 

ICG  10:1  Uptake  Ratio 
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Ratio 
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Heterogeneity: 
Surroundings: 
100:1  Uptake 

ICG 

DTTCI 

Ratio 
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Imaging  of  fluorescent  yield  and  lifetime 
from  multiply  scattered  light  reemitted  from 
random  media 


D.  Y.  Paithankar,  A.  U.  Chen,  B.  W.  Pogue,  M.  S.  Patterson,  and  E.  M.  Sevick-Muraca 
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The  feasibility  of  employing  fluorescent  contrast  agents  to  perform  optical  imaging  in  tissues  and  other 
scattering  media  has  been  examined  through  computational  studies.  Fluorescence  lifetime  and  yield 
can  give  crucial  information  about  local  metabolite  concentrations  or  environmental  conditions  within 
tissues.  This  information  can  be  employed  toward  disease  detection,  diagnosis,  and  treatment  if  non- 
invasively  quantitated  from  reemitted  optical  signals.  However,  the  problem  of  inverse  image  recon¬ 
struction  of  fluorescence  yield  and  lifetime  is  complicated  because  of  the  highly  scattering  nature  of  the 
tissue.  Here  a  light  propagation  model  employing  the  diffusion  equation  is  used  to  account  for  the 
scattering  of  both  the  excitation  and  fluorescent  light.  Simulated  measurements  of  frequency-domain 
parameters  of  fluorescent  modulated  ac  amplitude  and  phase  lag  are  used  as  inputs  to  an  inverse 
image-reconstruction  algorithm,  which  employs  the  diffusion  model  to  predict  frequency-domain  mea¬ 
surements  resulting  from  a  modulated  input  at  the  phantom  periphery.  In  the  inverse  image- 
reconstruction  algorithm,  a  Newton-Raphson  technique  combined  with  a  Marquardt  algorithm  is 
employed  to  converge  on  the  fluorescent  properties  within  the  medium.  The  successful  reconstruction 
of  both  the  fluorescence  yield  and  lifetime  in  the  case  of  a  heterogeneous  fluorophore  distribution  within 
a  scattering  medium  has  been  demonstrated  without  a  priori  information  or  without  the  necessity  of 
obtaining  absence  images.  ©  1997  Optical  Society  of  America 

Key  words:  Fluorescence  imaging,  biomedical  optics,  image  reconstruction,  fluorescence  lifetime, 
fluorescence  yield. 


1.  Introduction 

Over  the  past  decade,  several  investigators  have  ex¬ 
plored  the  use  of  exogenous  fluorescent  dyes  as  con¬ 
trast  agents  to  differentiate  diseased  and  normal 
tissues  from  noninvasive  or  endoscopic  optical  mea¬ 
surements.  The  diagnosis  of  burn  depth  following 
Indocyanine  Green  dye  administration1  and  the  de¬ 
marcation  of  neoplastic  tissues  following  intravascu¬ 
lar  porphyrin  dye  administration2-4  are  possible 
because  of  their  leakage  from  vessels  corrupted  by 
insult  or  disease.  The  concomitant  increase  in  fluo- 
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rescence  from  the  dye  can  be  detected  at  the  tissue 
surface  and  permits  the  detection  of  disease.  With 
the  development  of  dyes  and  photodynamic  agents 
that  excite  and  reemit  in  the  near-infrared  red  wave¬ 
length  regime,  the  noninvasive  detection  of  diseased 
tissues  located  deep  within  tissues  may  also  be  pos¬ 
sible  because  red  excitation  and  reemission  light  can 
travel  significant  distances  to  and  from  the  tissue-air 
interface.5  However,  a  long-standing  problem  has 
been  the  low  uptake  or  leakage  into  neoplastic  tis¬ 
sues,  providing  insufficient  contrast  for  the  detection 
of  diseased  tissues.  Although  targeted  delivery  may 
improve  uptake  ratios  and  contrast  imparted  by 
these  otherwise  therapeutic  agents,  these  approaches 
have  been  elusive  even  in  conventional  imaging  mo¬ 
dalities  such  as  magnetic  resonance  imaging  (MRI) 
and  x-ray  computed  tomography.  Investigators 
have  sought  to  alleviate  the  optical  contrast  uptake 
problem  by  employing  dyes  that  fluoresce  differently 
in  diseased  tissues  than  in  normal  tissues.  The  use 
of  agents  whose  reemission  characteristics  vary  with 
tissue  pH6’7  and  p028  may  not  only  provide  detection 
of  diseased  tissues  by  the  nature  of  differing  fluores- 
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cent  properties  but  may  also  contain  diagnostic  infor¬ 
mation.  The  difficulty  lies  in  measuring  the 
multiply  scattered  reemitted  light  and  reconstructing 
an  image  that  diffei'entiates  tissues  on  the  basis  of 
fluorescent  properties,  such  as  fluorescent  yield  and 
lifetime. 

In  this  study  we  investigate  the  combination  of 
photon  migration  imaging  techniques  with  fluores¬ 
cence  spectroscopy  to  examine  the  feasibility  of  fluo¬ 
rescence  lifetime  imaging  (FLI).9  FLI  may  be 
considered  an  optical  analog  to  MRI.  Whereas  MRI 
depends  on  monitoring  the  spatial  variation  in  the 
relaxation  time  of  spin-spin  states  to  provide  high- 
resolution  imaging  in  tissues,  FLI  depends  on  moni¬ 
toring  the  spatial  variation  in  the  fluorophore  lifetime 
and  yield.  The  complication  with  FLI  in  tissues 
arises  from  the  excitation  and  fluorescent  photon 
times  of  flight,  which  are  similar  in  magnitude  to  the 
fluorescence  lifetime.  Herein  we  combine  inverse 
imaging  techniques  (as  developed  previously  by 
Pogue  et  al.10)  and  fluorescence  lifetime  spectroscopy 
to  demonstrate  the  feasibility  of  FLI  by  using  com¬ 
putational  studies.  We  separately  introduce  the 
concepts  of  photon  migration  imaging  and  fluores¬ 
cence  spectroscopy  in  Section  2  and  then  describe  the 
combination  of  techniques  in  our  computational  ap¬ 
proach  to  FLI.  Our  reconstructed  images  show  the 
ability  to  image  both  fluorescence  lifetime  and  fluo¬ 
rescence  yield  independent  of  any  a  priori  informa¬ 
tion  and  point  to  the  development  of  fluorescent  dyes 
whose  reemission  characteristics  are  environmen¬ 
tally  sensitive  and  provide  contrast  for  the  optical 
detection  of  diseased  tissues. 

2.  Background:  Photon  Migration  Imaging  and 
Fluorescence  Spectroscopy 

The  ability  to  reconstruct  an  internal  map  of  absorp¬ 
tion  optical  properties  from  continuous  wave11  and 
absorption  and  scattering  properties  from  frequency- 
domain10’12-13  measurements  has  been  successfully 
demonstrated  by  using  multiply  scattered  light  from 
random  media.  Continuous  wave  measurements 
consist  of  noninvasively  monitoring  the  attenuation 
of  light  as  a  function  of  position  around  the  periphery 
of  a  heterogeneous  tissue  phantom  model  in  response 
to  a  constant  intensity  source,  whereas  frequency- 
domain  measurements  consist  of  monitoring  the 
phase  and  amplitude  modulation  of  multiply  scat¬ 
tered  light  at  various  peripheral  positions  in  response 
to  an  incident  modulated  light  source.  Both  ap¬ 
proaches  employ  a  perturbation  analysis  whereby  the 
absorption  and  scattering  properties  at  each  pixel 
position  within  the  tissue  phantom  are  individually 
and  independently  adjusted  until  the  measurements 
of  reemitted  light  at  the  periphery  of  the  tissue  phan¬ 
tom  match  those  predicted  by  a  model  for  diffusive 
light  propagation  in  random  media.  These  studies 
show  the  potential  for  image  reconstruction  of  hidden 
tissue  heterogeneities  when  the  contrast  caused  by 
either  absorption  or  scattering  coefficients  is  as  low  as 
2:1. 13  However,  Troy  et  al.14  conducted  in  vitro  mea¬ 
surements  of  normal  and  diseased  breast  tissues  and 


showed  that  the  endogenous  optical  contrast  between 
120  normal  and  diseased  breast  tissues  is  not  signif¬ 
icantly  different  for  consistent  detection  with  optical 
techniques.  These  measurements  were  conducted 
ex  vivo  and  consequently  may  underestimate  the  con¬ 
trast  available  in  vivo  caused  by  increased  tumor 
vasculature  and  absorption  caused  by  hemoglobin. 
Nonetheless,  these  results  suggest  the  need  for  opti¬ 
cal  contrast  agents. 

Previously,  we  experimentally  demonstrated  the 
increased  sensitivity  for  detecting  heterogeneities  on 
the  basis  of  fluorescence  as  opposed  to  absorption 
when  time-dependent  measurements  are  employed.15 
Indeed,  the  augmentation  of  optical  contrast  that  is 
due  to  the  lifetime  of  a  fluorophore  or  phosphorescent 
probe  has  been  recognized  by  Wu  et  al.16  in  tissue 
phantom  studies.  The  localization  of  a  fluorescing 
body  has  been  performed  by  various  researchers.17-19 
Upon  activation  into  a  higher  electronic  state  by  the 
absorption  of  light,  an  activated  fluorophore  may  un¬ 
dergo  nonradiative  decay  or  radiative  decay;  the  lat¬ 
ter  results  in  the  reemission  of  a  fluorescent  photon. 
The  yield  is  defined  by  the  fractional  number  of  flu¬ 
orescent  photons  reemitted  for  each  excitation  photon 
absorbed  or  the  fraction  of  decay  events  that  results 
in  the  emission  of  a  fluorescent  photon.  The  lifetime 
is  defined  as  the  mean  survival  time  of  the  activated 
fluorophore  or  the  mean  time  between  the  absorption 
of  an  excitation  photon  and  the  reemission  of  a  fluo¬ 
rescent  photon.  Because  the  stability  of  the  acti¬ 
vated  fluorophore  depends  on  its  environment,  both 
yield  and  lifetime  are  also  dependent  on  the  (bio) 
chemical  environment  in  which  the  fluorophore  re¬ 
sides.  Consequently,  the  lifetime  can  provide  con¬ 
trast  based  on  the  differences  in  biochemical 
environments  of  normal  and  diseased  tissues,  similar 
to  the  contrast  provided  by  relaxation  times  in  nu¬ 
clear  magnetic  resonance  imaging. 

The  use  of  lifetime  for  contrast  in  optical  imaging  in 
tissues  is  not  new.  In  studies  demonstrating  the 
noninvasive  differentiation  of  hematoporphyrin- 
laden  tumors  from  normal  tissue,  Cubeddu  et  al.20 
employed  the  comparatively  long  lifetime  (>15  ns)  of 
reemitted  fluorescence  caused  by  increased  porphy¬ 
rin  uptake  in  tumors  over  the  lifetime  of  endogenous 
compounds  (>7  ns)  in  normal  tissues  in  order  to  pro¬ 
vide  discrimination  of  the  two  tissue  types.  These 
microscopy  studies  involved  unscattered  light  that 
provided  correct  reemission  kinetics  of  measured  flu¬ 
orescence.  In  this  computational  study,  we  employ 
contrast  that  is  due  to  increased  uptake  as  well  as 
lifetime  in  order  to  differentiate  diseased  tissue  by 
using  multiply  scattered  light  reemitted  from  tissues 
that  involve  deeply  seated  tissue  volumes.  The  re¬ 
emission  kinetics  of  phosphorescent  probes  do  not 
permit  an  interrogation  deep  into  tissue  with  time- 
domain  approaches.21  Our  preliminary  frequency- 
domain  computations  (unpublished)  suggest  that  the 
simultaneous  determination  of  location  and  lifetime 
is  difficult  when  long-lived  phosphorescent  probes  are 
used.  In  Section  3  we  underscore  the  magnitude  of 
the  problem  by  describing  the  forward  imaging  prob- 
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lem  (i.e.,  prediction  of  the  frequency-domain  mea¬ 
surements  given  the  location,  optical  properties,  and 
fluorescent  properties  of  the  random  media)  as  well 
as  our  adapted  approach  for  solving  the  inverse  prob¬ 
lem  (i.e.,  prediction  of  the  location,  optical  properties, 
and  fluorescent  properties  of  the  medium  from  mea¬ 
surements  of  frequency-domain  light  propagation). 

3.  Forward  Problem  and  the  Solution  Methodology 

The  spatial  and  temporal  transport  of  light  in  tissues 
or  multiply  scattering  media  can  be  accurately  de¬ 
scribed  by  the  diffusion  approximation  to  the  radiative 
transport  equation.  A  coupled  frequency-domain 
diffusion  equation  can  be  used  to  predict  the  excita¬ 
tion  and  emission  fluence  rates,  <bx{r,  to)  and  4>m(r,  co), 
at  any  location  r  within  a  tissue  phantom  by  Eqs.  (1) 
and  (2),  respectively.22'24 

V  •  [Dx(r)V<Z>x(r,  w)]  -  [gjr)  +  i<a/cn] 

X  4 \{r,  oj)  +  Sx(r,  «)  =  0,  (1) 

V  •  [Dm(r)Vd)m(r,  w)]  -  [|x0m(r)  +  iw/cj 

X  4>m(r,  co)  +  Sm(r ,  co)  =  0.  (2) 

The  source  term  for  excitation  light  Sx(r,  co)  is  due  to 
the  sinusoidally  modulated  light  at  frequency  co  = 
2-it f,  where  f  is  usually  in  the  megahertz  range.  The 
first  term  in  both  equations  represents  the  diffusive 
or  random-walk  transport  of  light  where  Dxm  is  the 
optical  diffusion  coefficient,  i.e., 

Ae,m  =  [3(|J-QIm  +  As**/)]  S  (3) 

and  and  gs'  are  the  absorption  and  isotropic  scat¬ 
tering  coefficients,  respectively.  The  optical  proper¬ 
ties  are  dependent  on  the  wavelength  of  light  and 
thus  are  different  for  the  excitation  (subscript  x)  and 
fluorescent  (subscript  m)  light.  The  total  absorption 
coefficient  at  the  excitation  wavelength,  p,a  ,  is  due  to 
contributions  from  nonfluorescing  chromophores  as 
well  as  from  fluorescent  dye.  The  total  absorption 
coefficient  is  given  by  the  sum  of  absorption  coeffi¬ 
cients  that  are  due  to  nonfluorescing  chromophores, 
|xa  and  fluorophores,  \\.„  .  We  assume  that  the 

absorption  experienced  at  the  fluorescent  wavelength 
primarily  is  due  to  nonfluorescing  chromophores. 
The  velocity  of  light  in  tissue  is  cn  =  c/n,  where  n  is 
the  average  index  of  refraction.  The  source  term  for 
the  fluorescent  light  is  dependent  on  the  excitation 
light  fluence,  4>x(r,  co),  and  is  given  by 

1  —  tcoT(r) 

Sm(r,  co)  =  Tpr  (r)d)x(r,  co)  — — y—.2 .  (4) 

x  1  +  co  T(r) 

This  term  arises  from  the  Fourier  transform  of  the 
single-exponential  fluorescence  decay  term  (1/ 
T)exp(-t/T)  in  the  time  domain  following  an  incident 
pulse  of  excitation  light  where  t  is  the  fluorophore 
lifetime.  Here  -q  is  the  quantum  yield  and  the  ab¬ 
sorption  coefficient,  p  a  is  the  product  of  the  extinc¬ 
tion  coefficient,  loge  1(5,  and  the  concentration  of  the 
fluorophore  in  the  ground  state.  For  the  purposes  of 


this  study,  the  combined  product,  Tipa  ,  is  termed 
the  fluorescent  yield  and  is  proportional  to  the  gen¬ 
erated  fluorescence  fluence.  Note  that  multiexpo¬ 
nential  time  decay  can  also  be  handled  with  this 
procedure  by  a  simple  extension. 

In  the  source  term  for  fluorescent  light  given  in  Eq. 
(4),  changes  as  the  relative  fractions  of  the  flu¬ 
orophore  in  the  ground  and  excited  states  change. 
The  saturation  effects  would  have  to  be  handled  by 
taking  into  account  the  relative  amounts  of  fluoro¬ 
phore  in  the  ground  and  excited  states.  We  neglect 
the  saturation  effects  and  assume  single-exponential 
decay  kinetics  in  this  initial  study. 

Furthermore,  pa  in  Eq.  (1)  is  the  sum  of  |xti  and 
pn  .  In  what  follows,  images  of  lifetime,  t,  and 
images  of  Tip ,  are  obtained.  There  is  a  difficulty 
in  the  estimation  of  p^  that  is  given  by  p.„  =  p  + 
pa  _  ,  because  the  explicit  values  of  pa  are  not 
known  (the  values  of  the  product  of  t)  and  p„  are 
known).  Because  the  major  contribution  to  ga  is 
from  pa  and  not  pa  _  ,  we  have  used  an  approxi¬ 
mate  expression  for  pa  ,  pa^  =  \\,u  +  T)pa  m  and  have 

chosen  r]  as  1  only  for  the  purposes  of  calculation  of 
pa  .  Our  approach  can  be  extended  to  permit  imag¬ 
ing  of  p,  pa  ,  pa  _,  and  t  from  measurements  con¬ 
ducted  at*  botfi  excitation  and  fluorescent 
wavelengths.  Although  reconstruction  from  fluores¬ 
cent  measurements  provides  r)pa  __  ,  excitation  mea¬ 
surements  can  provide  pa  ^  (in  tne  absence  of  a 
fluorophore)  and  pa  _  +  pa  (in  the  presence  of  a 
fluorophore).  From  these  images,  maps  of  t]  and 
pQ  can  be  obtained  in  principle. 

Both  Eqs.  (1)  and  (2)  are  linear  complex  elliptic 
equations  that  can  be  solved  as  boundary  value 
problems  for  complex  quantities  <&x(r,  co)  and  4>m(r, 
co).  We  employ  the  method  of  finite  differences  in 
which  we  place  a  grid  over  the  space  domain  and 
obtain  an  approximation  to  the  solution  at  each  grid 
point,  j.  One  of  the  fastest  methods  to  solve  these 
linear  elliptic  boundary  value  problems  is  the  mul¬ 
tigrid  solution  (see  the  review  by  Fulton  et  al.25). 
In  the  procedure  for  the  multigrid  solution,  an  ini¬ 
tial  solution  is  obtained  quickly  for  coarse  grids  that 
is  then  further  refined  for  a  better  solution  for  finer 
grids.  This  is  an  involved  procedure  in  which  we 
have  elected  to  use  mudpack  routines.26  mudpack 
routines  are  flexible  and  allow  placement  of  sources 
either  at  the  surface  or  inside  the  phantom.  For 
the  equations  to  be  solved,  it  is  assumed  that 
w)  =  0  on  the  tissue  surface,  which  is  known  as  the 
zero-fluence  boundary  condition.  This  is  imple¬ 
mented  by  assigning  the  absorption  coefficient  for 
both  excitation  and  fluorescent  light  at  all  the  grid 
points  in  the  square  grid  lying  outside  the  circular 
tissue  phantom  to  a  large  value.  The  source  was 
simulated  by  setting  the  value  of  <1>X  to  an  arbitrary 
complex  number  at  a  grid  point  on  the  surface 
where  the  source  is  located.  The  solution  of  Eqs. 
(1)  and  (2)  yields  a  complex  number  for  4>m  at  each 
grid  point,  j.  The  detected  signal  at  the  surface  is 
proportional  to  the  normal  component  of  the  gradi- 
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Fig.  1.  Schematic  of  the  circular  simulated  tissue  phantom  inter¬ 
rogated  by  four  sources,  one  source  at  a  time.  Twenty  detectors 
are  located  on  the  periphery  with  uniform  separation.  A  circular 
object  that  mimics  a  hidden  diseased  tissue  volume,  located  within 
the  tissue  phantom,  is  also  shown. 


ent.  In  this  study,  for  evaluation  of  the  signal  at 
detector  i  located  on  the  surface,  we  used  the  value 
of  <hm  at  an  internal  grid  point  closest  to  the  detec¬ 
tor.  This  is  reasonable  in  our  initial  approach  be¬ 
cause  the  normal  component  of  the  gradient  is 
proportional  to  <£>m  just  inside  the  surface.10  The 
phase  lag,  0M,  and  the  log  of  the  ac  amplitude,  Mm, 
at  the  detectors  were  calculated  with  respect  to  the 
phase  and  the  ac  amplitude  of  the  source.  The 
computational  time  required  to  solve  for  the  fluo¬ 
rescent  fluence  on  a  SunSparc  10  for  a  17  X  17  grid 
is  0.04  s,  and  that  for  a  33  X  33  grid  is  0.165  s. 
Different  grid  sizes  were  used  in  the  testing  phase  of 
the  solution  of  the  equations,  and  the  solutions  were 
in  agreement  within  the  numerical  error  caused  by 
the  finite  grid  size. 

Before  attempting  to  solve  the  inverse  problem  for 
the  simulated  phantom  illustrated  in  Fig.  1,  we  must 
first  understand  the  effects  of  changing  the  fluores¬ 
cent  optical  properties  of  the  tissue  on  0m  and  Mm 
measured  at  a  detector  or  series  of  detectors  located 
on  its  surface.  Solutions  to  Eqs.  (1)  and  (2)  were 
obtained  in  two  dimensions  for  a  65  X  65  grid  cover¬ 
ing  a  100-mm-diameter  circular  tissue  phantom  with 
a  circular  embedded  heterogeneity  of  30  mm  diame¬ 
ter  and  located  at  the  center  of  the  tissue  phantom. 
The  simulated  measurements  of  fluorescent  phase 
shift  and  ac  amplitude  are  reported  for  20,  equally 
spaced,  circumferentially  located  detectors.  The 
modulation  frequency,  f,  was  set  equal  to  150  MHz. 
The  optical  properties  of  the  heterogeneity  and  the 
background  are  shown  in  Table  1. 


Fig.  2.  Log  plot  of  fluorescent  light  ac  amplitude  at  150  MHz  at 
various  detectors  for  a  heterogeneous  tissue  phantom  with  a  30- 
mm-diameter  object  located  at  the  center  of  the  phantom  for  var¬ 
ious  object  np.a  values.  Absorption  coefficient  up.,,  is  1  X 
10-5  mm-1  for  the  background;  lifetimes  t  for  both  the  object  and 
the  background  are  1  ns. 


A.  Effect  of  Tn|J-ax_,m  0n  and  Mm 
In  order  to  evaluate  the  influence  of  r)|jLa  ,  0m  and 
Mm  were  computed  at  each  detector  as  the  value  of 
T)[jL0t  in  the  heterogeneity  increased  from  10  4 
mrrT  1  to  lO-1  mm-1  and  as  ri(xa  in  the  background 
was  maintained  constant.  The  modulation  fre¬ 
quency  was  chosen  as  150  MHz.  The  lifetime,  t,  was 
set  equal  to  1  ns  for  both  the  object  and  the  back¬ 
ground  causing  contrast  resulting  solely  from  differ¬ 
ences  in  T||xaj  .  The  plots  of  0m  and  Mm  are  shown 
in  Figs.  2  and  3,  respectively,  for  one  source  located  at 
position  A.  The  curves  for  Mm  are  not  smooth  be¬ 
cause  the  circular  surface  is  approximated  by  joining 
discrete  grid  points  closest  to  the  circumference  by 
straight  lines.  Furthermore,  the  curves  for  both  0m 
and  Mm  are  symmetric  for  detectors  2  and  20,  3  and 
19,  4  and  18,  and  so  on  because  these  detector  pairs 
are  symmetric  with  respect  to  the  source.  The  de¬ 
tectors  near  the  source  are  unaffected  by  the  presence 
of  an  object  because  the  photon  sampling  volume  does 
not  include  the  region  occupied  by  the  heterogeneity. 
At  the  other  detector  locations,  the  fluorescent  ac 
amplitude  increases  as  r\\xa^m  of  the  heterogeneity 
increases.  As  _  of  the  heterogeneity  increases 
to  high  values,  the  ac  amplitude  approaches  an  upper 
limit  that  is  due  to  high  absorption  coefficient  ga  , 
which  is  a  sum  of  p,  and  g0  .  This  is  similar  to 
the  inner-filter  effect  in  which  the  high  absorption  in 
the  heterogeneity  shields  the  interior  of  the  hetero¬ 
geneity  from  excitation  photons.27  Figure  3  illus¬ 
trates  changes  in  fluorescent  phase  shift,  0m,  as  a 


Table  1.  Optical  Properties  and  Experimental  Parameters  for  the  Forward  Problem 


Fa, 

A,,'  or  |xSm' 

Fa,_ 

•niJ-a,.,,  Background 

t  Background 

Frequency 

(mm-1) 

(mm"1) 

(mm-1) 

(mm-1) 

(mm-1) 

(ns) 

(MHz) 

0.0 

1.0 

0.0 

1.0  X  10-5 

1.0 

150.0 
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Fig.  3.  Plot  of  fluorescent  light  phase  shift  at  150  MHz  at  various 
detectors  for  a  heterogeneous  tissue  phantom  with  a  30-mm- 
diameter  object  located  at  the  center  of  the  phantom  for  various 
object  T||i0  _  values.  Absorption  coefficients  and  lifetimes  are  the 
same  as  in  Fig.  2. 


E 

< 


Fig.  4.  Log  plot  of  fluorescent  light  ac  amplitude  at  150  MHz  at 
various  detectors  for  a  heterogeneous  tissue  phantom  with  a  30- 
mm-diameter  object  located  at  the  center  of  the  phantom  for  var¬ 
ious  object  t  values.  Background  lifetime  t  is  1  ns  and  absorption 
coefficients  T||ia  ^  are  1  x  10-5  mnri1  for  the  background  and  1  X 
10-3  mm-1  for  the  object. 


function  of  the  fluorescent  yield  that  is  due  to  the 
fluorophore,  -ijp.^  As  T||xa  _  is  increased  from 
10-3  mm-1  to  10-1  mm-1  (when  the  background 
value  is  10-5  mm-1),  little  change  in  phase  shift  oc¬ 
curs.  Because  the  phase  shift  of  reemitted  fluores¬ 
cence  in  dilute,  nonscattering  solutions  is  a  function 
of  lifetime  and  is  not  dependent  on  concentration  [0  ~ 
tan-1(coij],  one  might  expect  similar  trends  in  scat¬ 
tering  media.  Indeed,  as  shown  in  Fig.  3,  phase- 
shift  changes  caused  by  fluorophore  concentration 
differences  are  small.  Consequently,  the  changes  in 
phase  shift  with  increasing  irpq,  m  of  the  heterogene¬ 
ity  can  be  attributed  to  the  alteration  in  photon  mi¬ 
gration  as  a  result  of  the  presence  of  heterogeneities 
with  a  high  absorption  of  excitation  light.28  The 
presence  of  heterogeneities  with  high  absorption  re¬ 
duces  the  effective  path  length  of  photon  migration 
and  reduces  the  phase  shift,  though  this  reduction  is 
small  in  magnitude.  In  summary,  Mm  is  directly 
and  strongly  dependent  on  changes  in  rpx,,  ^  of  a 
simulated  tissue  heterogeneity,  whereas  0m  is  indi¬ 
rectly  and  weakly  dependent  on  r||xa  m  because  of 
changes  in  photon  migration. 

B.  Effect  of  t  on  Mm  and  6m 

In  order  to  evaluate  the  influence  of  t,  Mm  and  0m 
were  calculated  at  each  detector  as  the  values  of  t  in 
the  heterogeneity  varied  from  10-1  ns  to  103  ns  and 
the  value  of  t  in  the  background  was  held  at  1  ns. 
The  modulation  frequency  was  at  150  MHz.  Back¬ 
ground  T|(jia^  m  was  set  to  10-5  mm-1,  and  r\\i.a  for 
the  object  was  set  to  10-3  mm-1.  As  shown  in  Fig.  4, 
the  detected  ac  amplitude  increases  as  t  decreases. 
Because  Mm  of  reemitted  fluorescence  in  a  dilute, 
nonscattering  solution  is  a  function  of  lifetime  (and 
also  fluorophore  concentration),  one  would  expect 
that  in  the  presence  of  scatter,  similar  trends  would 
be  observed.  Figure  5  illustrates  the  values  of  the 
fluorescent  phase  shift  at  each  detector  as  the  life¬ 
time  of  the  heterogeneity  is  changed  from  0.1  ns  to 
1000  ns.  At  the  given  modulation  frequency  (150 


MHz  in  this  calculation),  Qm  first  increases,  reaches  a 
maximum,  and  then  subsequently  decreases  as  t  is 
increased  from  0.1  ns  to  1000  ns.  There  is  a  complex 
interplay  of  the  signal  coming  from  the  object  for 
which  we  are  changing  the  lifetime  and  the  signal 
from  the  background  that  leads  to  the  behavior  de¬ 
scribed.  In  summary,  both  Mm  and  0m  at  each  de¬ 
tector  are  directly  and  strongly  influenced  by  the 
value  of  lifetime  in  the  heterogeneity. 

Detailed  analytical  expressions  for  fluorescence 
amplitude  and  phase  shift  for  infinite  and  semi¬ 
infinite  media  with  spherical  heterogeneities  have 
been  provided  by  Li  et  al.29  Our  numerical  compu¬ 
tational  methods  provide  a  simulation  of  finite  media 
and  arbitrary  shaped  hidden  objects.  Our  numeri¬ 
cal  results  agree  with  the  general  predictions  pro¬ 
vided  by  Li  et  al.29 

The  solution  of  the  forward  problem  was  used  as 
inputs  to  the  solution  of  the  inverse  problem  de- 


Fig.  5.  Plot  of  fluorescent  light  phase  shift  at  150  MHz  at  various 
detectors  for  a  heterogeneous  tissue  phantom  with  a  30-mm- 
diameter  object  located  at  the  center  of  the  phantom  for  various 
object  t  values.  Background  lifetime  t  is  1  ns  and  absorption 
coefficients  Tip.,^  m  are  1  X  10-5  mm-1  for  the  background  and  1  X 
10“3  mnT1  for  the  object. 
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Table  2.  Optical  Properties  Used  to  Generate  the  Simulated  Experimental  Data  as  Inputs  to  the  Inverse  Image-Reconstruction  Algorithm 


Case 

(mm-1) 

Pam 

(mm-1) 

(mm-1) 

T 

Background 

(ns) 

T)|Xa 

1 1  c ix—*m 

Background 

(mm-1) 

uM  Gaussian  Noise 
in  Log  of  ac 
Amplitude 

ae  Gaussian  Noise 
in  Phase 
(deg) 

5.A 

0.0 

0.0 

1.0 

10.0 

1.0  X  10-5 

0.01 

0.1 

5.B 

1.0  X  10-3 

1.0  X  10-3 

1.0 

10.0 

1.0  X  10-5 

0.01 

0.1 

5.C 

0.0 

0.0 

1.0 

10.0 

1.0  X  10-5 

0.01 

1.0 

5.D 

1.0  x  10-3 

1.0  X  10-3 

1.0 

10.0 

1.0  X  10-5 

0.01 

0.1 

scribed  in  Section  4.  Forward  solutions  in  two  di¬ 
mensions  were  computed  for  the  three  cases  outlined 
in  Tables  2  and  3  for  the  source  and  detector  geom¬ 
etry  in  a  100-mm-diameter  simulated  tissue  phantom 
(Fig.  1).  Values  of  0m  and  Mm  were  computed  at 
each  detector  in  response  to  four  sources  of  excitation 
light  located  at  the  periphery  modulated  at  a  single 
frequency  of  150  MHz.  Consequently,  the  forward 
solution  predicted  80  values  of  9m  and  Mm,  each  at 
detector  i  (i  =  1,  20)  in  response  to  source  k  (k  =  1, 4). 
Gaussian  noise  with  a  standard  deviation  of  0.1°  (or 
a  liberal  1°)  in  0m  and  1%  in  Mm  was  superimposed  on 
the  solution  of  the  forward  problem.  These  obtained 
data  sets  are  used  as  the  simulated  experiments  as 
input  to  the  reconstruction  algorithm  described  be¬ 
low. 


ulated  experimental  values  of  (0m),  and  (Mm)f  re¬ 
sulting  from  each  source,  k,  k  =  1,  4.  The  choice  of 
appropriate  functions  to  be  minimized  for  our  re¬ 
constructions  is  now  discussed.  Section  3  de¬ 
scribes  at  length  how  changes  in  the  fluorescent 
yield  and  lifetime  affect  the  measured  values  of  ac 
magnitude  and  phase.  Changes  in  fluorescent 
yield,  Tpx0  ,  have  a  direct  and  strong  influence  on 
the  ac  magnitude  and  only  an  indirect  and  weak 
influence  on  the  phase.  Hence  we  have  chosen  to 
adjust  iteratively  the  values  of  ("mv™);  in  order  to 
minimize  the  merit  function,  x,x2>  that  depends  on 
the  ac  magnitude  measurements: 


X 


2 

M- 


i  4  i  20  Im  -m  Y 
4&20&  \  <yM  ) 


(5) 


4.  Inverse  Problem  and  the  Solution  Methodology 

We  begin  our  inverse  calculations  from  a  uniform 
starting  guess  of  a  fluorescence  yield,  (imxa  )p  and 
lifetime,  (t),-,  at  all  the  grid  points  within  the  tissue 
phantom.  With  the  initial  guess  of  an  optical  prop¬ 
erty  map,  one  can  evaluate  the  complex  fluence  at 
the  various  detector  locations.  There  are  two  ap¬ 
proaches  that  one  can  follow:  one  involves  mini¬ 
mizing  the  difference  between  the  complex  number 
that  is  calculated  from  the  ac  magnitude  and  phase, 
and  the  other  uses  the  ac  magnitude  and  phase 
information  itself.  The  two  approaches  are  equiv¬ 
alent,  and  we  have  chosen  the  second  approach  of 
using  the  ac  magnitude  and  phase.  With  the  start¬ 
ing  guess  of  optical  properties,  we  compute  a  pre¬ 
diction  of  phase  shift  (0m)j  and  log  of  ac  amplitude 
(Mm)i  at  each  detector  i.  Values  of  (q|jia  )•  and 
(t)7-  are  iteratively  adjusted  at  each  grid  point  to 


where  txM  is  the  typical  standard  deviation  of  noise  in 
Mm,  taken  to  be  0.01. 

The  value  of  t  is  the  simulated  experimental 
value  computed  in  Section  3  with  added  Gaussian 
noise.  The  goal  of  the  algorithm  is  to  minimize  x(i2 
by  appropriate  updates  of  (r)|xa After  an  itera¬ 
tion  of  the  (irip.a  )j-  update,  values  of  (t )•  were  ad¬ 
justed  in  the  next  iteration  in  order  to  minimize  a 
second  merit  function,  \r2: 


-l  4  -1  20 

Xt  =4j§2()|' 


^mobs,i  +  (®mobai  ~  0m, i 


% 


UR 


(6) 


where  ue  is  the  typical  standard  deviation  of  noise 
in  0,  taken  to  be  1°.  The  values  of  0m  b=  are  the 
simulated  experimental  values  computed  in  Section 


minimize 

the  error  between  the  predicted  and  sim-  3 

Table  3.  Location  and  Area  of  the  Simulated  Heterogeneities: 

with  added  Gaussian  noise.  The  above  merit 

Comparison  of  Expected  and  Reconstructed  Values 

Object  1 

Object  2 

Area 

Location  (x,  y) 

Area 

Location  ( x ,  y) 

Case 

(mm2) 

(mm) 

(mm2) 

(mm) 

5.A 

706.0  (expected) 

(60,  60)  (expected) 

not  applicable 

not  applicable 

742  (obtained) 

(61,  59)  (obtained) 

5.B 

706.0  (expected) 

(60,  60)  (expected) 

not  applicable 

not  applicable 

683  (obtained) 

(59,  58)  (obtained) 

5.C 

314.1  (expected) 

(32.3,  67.7)  (expected) 

314.1  (expected) 

(67.7,  32.3)  (expected) 

381  (obtained) 

(34,  68)  (obtained) 

342  (obtained) 

(65,  35)  (obtained) 

5.D 

706.0  (expected) 

(60,  60)  (expected) 

not  applicable 

not  applicable 

693  (obtained) 

(61,  57)  (obtained) 
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function  is  chosen  because  changes  in  lifetime  t 
have  a  direct  and  strong  influence  on  both  the  ac 
magnitude  and  phase.  Hence  we  have  chosen  to 
adjust  iteratively  the  values  of  (t);-  in  order  to  min¬ 
imize  the  merit  function  that  consists  of  both  ac 
magnitude  and  phase  measurements.  The  goal  of 
this  part  of  the  algorithm  is  to  minimize  yT2  by 
appropriate  choice  of  (t)j  as  described  below.  In 
Eq.  (5),  the  choice  of  the  value  of  uM  does  not  affect 
the  minimization- update  process.  We  have  cho¬ 
sen  the  typical  noise  values  for  M  and  0,  crM  and  o-9, 
respectively,  to  be  factors  for  converting  to  dimen¬ 
sionless  numbers.  In  fact  uM  and  ae  are  simply 
scale  factors  that  assess  confidence  in  measure¬ 
ment.  For  example,  when  ct6  is  smaller  for  a  given 
ctm,  the  0  residuals  are  weighted  more  heavily  in  the 
inversion  as  compared  with  the  M  residuals. 
Thus,  factor  a  is  necessary  to  normalize  the  indi¬ 
vidual  sum  of  residuals  of  M  and  0  so  the  two  can  be 
added  in  Eq.  (6). 

To  update  values  of  (p|xa  X-  and  (t)j,  one  needs  the 
Jacobian  matrices  that  describe  the  sensitivity  of  the 
detector  response  at  position  i  to  changes  in  (t),-  and 
(p  each  grid  point,  j.  The  elements  of  the 

three  Jacobian=  matrices  employed,  J(M,  p(xa  ), 
J(M,  t),  and  J(0,  t),  are  given  by  jLj  =  [dMJ/ 
=  (dMJdTj),  and jtJ  =  {dQjdij),  respec¬ 
tively.  One  calculates  these  elements  by  solving  the 
forward  problem  four  times  for  each  grid  point,  j  to 
obtain  Mm  i  and  0m  i-  calculated  with  (t)7-  and  (t  +  At); 
and  with  (pqa^m)7'  and  (pix^  ;  +  Aq|r„rJ;.  From 
the  least-squares  minimization,  one  can  show  that 
the  update  in  and  t  can  be  calculated  by  using 

an  algorithm  (Newton’s  method)  similar  to  that  sug¬ 
gested  by  Yorkey  et  al.30  for  the  reconstruction  of 
images  obtained  by  electrical  impedance  tomogra¬ 
phy. 

Equations  (7)  and  (8)  provide  updates 


J(M,  pip,  pp. 


°Af 


+  Xj/ 


j{M,  pp0r_.J7 


°Jkf 


[Ap^OMm] 


,  (7) 


J(M,  T f j(M,  T)  ,  j(0,  T)rJ(0,  T)  ,  ^  = 
9  +  2  '  ^2* 


j(M,  T  f 


[At] 


J(0,  t)t  - 

2  W-mohs  ~  Mm)  H  2  (®«Jobs  —  9m) 


0"e 


(8) 


[Apiv^J,  and  [At]  to  the  yield  [p]Xa  _  ]  and  lifetime 
[t]  vectors  at  each  iteration.  The  update  of  [ppQ„„J 
is  based  on  the  ac  amplitude  data,  whereas  that  of  [t] 
is  based  on  both  the  ac  amplitude  and  phase  data. 
This  follows  from  the  discussion  in  the  previous  sec¬ 
tion  on  the  forward  problem.  and  Mm  are  the 

experimentally  simulated  and  calculated  vectors 
consisting  of  the  log  of  the  ac  amplitude  at  each 
of  the  i  detectors,  respectively.  Because  of  the  ill- 


Fig.  6.  Graphs  depicting  the  convergence  of  (a)  T|p.ai_m,  (b)  t  ver¬ 
sus  the  number  of  iterations  during  the  image  reconstruction  for 
the  case  described  in  Subsection  5.A.  Convergence  is  seen  to  be 
achieved  within  20  iterations  for  (a)  and  50  iterations  for  (b). 


conditioned  nature  of  the  Jacobian  matrices,  terms 
\jl  or  \2I  are  added  ( I  is  an  identity  matrix)  to  make 
the  matrices  more  diagonally  dominant  and  the  so¬ 
lution  of  the  algebraic  equations  more  robust.  Pa¬ 
rameters  X3  or  \2  are  adjusted  by  means  of  a 
Marquardt-Levenberg  type  of  algorithm.31  Lower 
triangular- upper  triangular  decomposition  and 
backsubstitution  are  employed  to  solve  the  simulta¬ 
neous  linear  algebraic  equations  given  in  Eqs.  (7)  and 
(8).31  At  each  iteration,  the  merit  functions,  the  Ja¬ 
cobian  matrices,  and  the  updates  in  the  fluorescent 
yield  and  lifetime  are  evaluated  and  the  iterations 
are  continued  until  the  convergence  criterion  is  met. 
Convergence  is  achieved  when  any  of  the  following 
three  quantities,  i.e.,  y2,  change  in  y2  in  successive 
iterations,  and  relative  change  in  y2  in  successive 
iterations,  is  lower  than  a  predetermined  value  of 
1.0  X  10“5. 

5.  Results  and  Discussion 

The  performance  of  FLI  by  using  the  inversion  algo¬ 
rithms  described  above  is  depicted  in  Figs.  6-9  for  the 
case  studies  listed  in  Tables  2  and  3.  Simulated 
experiment  5. A  was  designed  to  reconstruct  (t)7  and 
(pip,  Jj  with  no  absorption  that  is  due  to  nonfluo- 
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x  coordinate  (mm) 

(a) 


(b) 

Fig.  7.  Reconstructed  spatial  map  of  fluorescence  (a)  yield, 
AM Ob)  lifetime,  t,  on  a  two-dimensional,  17  X  17  grid  for  the 
case  described  in  Subsection  5.A.  |jvl  at  the  excitation  wave¬ 
length  not  accounting  for  fluorophore  absorption  is  zero;  \x.am  at  the 
emission  wavelength  is  also  zero.  The  average  values  of  np-a 
and  t  within  the  object  were  1  X  10~3  mirO1  and  1  ns,  respectively 
(expected)  and  0.93  X  10“3  mm-1  and  1.03,  respectively  (recon¬ 
structed).  Spurious  unphysically  high  values  of  np,0  and  t 
have  been  replaced  by  the  average  background  fluorescence  yield 
and  lifetime,  respectively,  obtained  from  the  inversion. 


rescing  chromophores,  whereas  simulated  experi¬ 
ment  5.B  included  a  finite  chromophore  absorption  of 
excitation  and  fluorescent  light  that  could  be  consid¬ 
ered  physiological.  Finally,  simulated  experiment 
5.C  evaluated  the  ability  to  determine  the  location  of 
two  hidden  objects  within  the  tissue  phantom, 
whereas  5.D  examined  the  reconstruction  when  the 
uptake  ratio  (ratio  of  fluorescent  yield  in  heterogene¬ 
ity  to  that  in  the  background)  was  20. 

A.  Single  Heterogeneity  with  Optical  Contrast  from  (t),- 
and  (TiM-a^  with  Absorption  of  Excitation  and  Fluorescent 
Light  by  Chromophores 

To  calculate  the  experimental  data  for  the  first  case, 
the  fluorescence  yields,  up,  ,  for  the  background 


x  coordinate  (mm) 

(a) 


x  coordinate  (mm) 

(b) 

Fig.  8.  Reconstructed  spatial  map  of  fluorescence  (a)  yield,  Tip^  , 
(b)  lifetime,  t,  on  a  two-dimensional,  17  X  17  grid  for  the  case 
described  in  Subsection  5.B.  |i„  at  the  excitation  wavelength  not 
accounting  for  fluorophore  absorption  is  1  X  1(T3  mm-1;  p.0m  at  the 
emission  wavelength  is  zero.  The  average  values  of  and  t 

within  the  object  were  1  X  10-3  mm-1  and  1  ns,  respectively  (ex¬ 
pected)  and  0.8  x  10-3  mm-1  and  0.7  ns,  respectively  (reconstruct¬ 
ed).  Spurious  unphysically  high  values  of  Tip,a  and  t  have  been 
replaced  by  the  average  background  fluorescence  yield  and  lifetime, 
respectively,  obtained  from  the  inversion. 

and  a  single  object  were  chosen  as  1  X  10”5  mm”1 
and  1  X  10”3  mm”1,  respectively,  and  the  lifetimes, 
t,  for  the  background  and  the  object  were  chosen  as 
10  and  1  ns,  respectively.  During  the  inverse  im¬ 
age  reconstruction,  no  a  priori  knowledge  of  either 
the  object  location  or  the  background  fluorescence 
properties  was  assumed  and  a  uniform  guess  of  1  X 
10”5  mm”1  and  10  ns  was  given  for  T)p,a  and  t, 
respectively.  Convergence  was  achieved  in  fewer 
than  50  iterations  (computational  time  on  a  Sun- 
Sparc  10  was  2  h)  for  a  two-dimensional  17  X  17 
grid.  The  average  values  of  r|p,a  and  t  in  the  grid 
points  that  occupy  the  simulated'  object  converge 
within  50  iterations  to  the  values  of  ^  =  0.93  X 
10”3  mm”1  and  t  =  1.03  ns,  which  are  close  to  the 
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Fig.  9.  Reconstructed  spatial  map  of  fluorescence  yield, 
on  a  two-dimensional,  33  x  33  grid  for  the  case  described  in  Sub¬ 
section  5.C.  at  the  excitation  wavelength  not  accounting  for 

fluorophore  absorption  is  zero;  p.0m  at  the  emission  wavelength  is 
also  zero.  The  Gaussian  noise  that  was  introduced  in  the  phase 
had  a  standard  deviation  of  1°.  The  object  locations  and  sizes  are 
recovered  correctly.  The  average  values  of  T)p.„  m  within  the  two 
objects  were  top  left,  1  X  10-3  mm-1  (expected)  and  2  X  10-3 
mm-1  (reconstructed);  bottom  right,  2  X  10-3  mm-1  (expected) 
and  1.8  X  10-3  mm-1  (reconstructed).  Spurious  unphysically 
high  values  of  T)p.a  ^m  have  been  replaced  by  the  average  back¬ 
ground  fluorescence  yield  obtained  from  the  inversion. 


correct  values  of  =  1  X  10*3  mm*1  and  t  = 

1  ns,  as  shown  in  Figs.  6(a)  and  6(b)  and  listed  in 
Table  4.  Figures  7(a)  and  7(b)  illustrate  the  recon¬ 
structed  images  of  and  t,  respectively,  and 

are  representative  of  the  expected  images.  The 
images  were  smoothed  by  interpolation  in  this  and 
subsequent  cases.  During  the  reconstruction,  it 
was  observed  that  some  grid  points,  often  on  or 
close  to  the  periphery,  had  unphysically  high  values 
of  the  yield  as  well  as  lifetime.  Typically,  these 
grid  points  were  loners  and  were  surrounded  by 
grid  points  with  reasonable  yield  and  lifetime  val¬ 
ues.  The  high  values  of  lifetime  lead  to  a  lower 


Table  4.  Fluorescence  Lifetime  t  and  Yield  m  of  the  Simulated 
Heterogeneities:  Comparison  of  Expected  and  Reconstructed  Values 


Case 

■HIV™,  (Object) 
(mm-1) 

t  (Object) 

(ns) 

5.A 

1.0  X  10-3  (expected) 

1.0  (expected) 

0.93  X  10-3  (obtained) 

1.03  (obtained) 

5.B 

1.0  X  10-3  (expected) 

1.0  (expected) 

2.1  X  10-3  (obtained) 

4.1  (obtained) 

5.C 

Top  left  object 

1.0  X  10-3  (expected) 

1.0  (expected) 

2  X  10-3  (obtained) 

4.1  (obtained) 

Bottom  right 

2.0  X  10-3  (expected) 

2.0  (expected) 

object 

1.8  X  10-3  (obtained) 

3.5  (obtained) 

5.D 

2.0  X  10-4  (expected) 

1.0  (expected) 

7.0  X  10-4  (obtained) 

4  (obtained) 

magnitude  of  fluorescence  generation  and  offset  the 
effect  of  high  yield  values.  We  believe  that  addi¬ 
tional  criteria,  such  as  smoothness  of  the  recon¬ 
structed  maps,  may  alleviate  the  spurious 
background  peak  values.  In  this  study,  these  spu¬ 
rious  values  were  replaced  by  the  average  back¬ 
ground  fluorescence  yield  and  lifetime  obtained 
from  the  inversion.  In  the  subsequent  reconstruc¬ 
tions,  spurious  values  were  similarly  handled. 

The  average  values  of  in  the  grid  points 

that  occupy  the  simulated  background  converge 
within  50  iterations  to  9  X  10“5  mm*1  in  compar¬ 
ison  with  the  correct  value  of  1  X  10*5  mm*1  (data 
not  shown  for  brevity).  The  value  of  background  t 
converges  to  5.4  ns,  which  is  also  far  from  the  cor¬ 
rect  value  of  10  ns.  Both  discrepancies  are  attrib¬ 
uted  to  the  fact  that  most  of  the  signal  contribution 
is  due  to  the  object,  with  the  object-to-background 
uptake  ratio  being  100:1.  The  dependence  of  the 
final  images  on  the  choice  of  the  initial  guess  was 
examined  by  providing  an  initial  uniform  guess  of 
1  X  10*4  mm-1  and  10  ns  for  T|ga  and  t,  respec¬ 
tively.  This  resulted  in  images  similar  to  those 
obtained  with  the  first  guess. 

The  location  of  the  heterogeneity  was  identified  as 
consisting  of  all  the  grid  points  with  t||jl„  higher 
than  35%  (arbitrarily  chosen)  of  the  peak  value  of 
T||xa  ^  [Fig.  7(a)].  The  average  of  the  coordinates  of 
all  the  identified  object  grid  points  was  the  position 
(60.8, 58.5),  which  is  close  to  position  (60, 60)  that  was 
used  to  simulate  the  experimental  data.  As  listed  in 
Table  3,  the  area  of  the  heterogeneity  based  on  our 
arbitrary  definition  for  identification  was  742  mm2, 
close  to  that  used  to  generate  our  simulated  experi¬ 
mental  data.  In  addition,  upon  inspection  of  Fig. 
7(a),  one  can  note  that  values  of  tqix„x  for  the  object 
on  grid  points  closer  to  the  center"  of  the  phan¬ 
tom  were  higher  than  those  toward  the  periphery. 
This  is  due  to  a  reduction  of  signal  contributions  of  grid 
points  located  farthest  away  from  source  and  detec¬ 
tors.  This  is  reflected  in  a  smaller  value  in  the 
Jacobian  matrices  and  results  in  a  poorer  reconstruc¬ 
tion  at  the  center  of  the  phantom. 

B.  Single  Heterogeneity  with  Optical  Contrast  from  (t)7- 
and  (t)  p.a  )/  with  Absorption  of  Excitation  and  Fluorescent 
Light  by  Chromophores 

The  same  hidden  object  location  as  well  as  optical 
parameters  were  used  as  described  in  Subsection 
5. A,  except  that  a  uniform  background  chro- 
mophore  absorption  coefficient  of  the  excitation 
light,  p,  of  1  X  10  3  mm  ,  and  of  the  fluorescent 
light,  pa'  of  1  X  10“3  mm*1,  was  used  to  generate 
the  simulated  experimental  data.  Although  exci¬ 
tation  light  propagation  was  not  employed  for  im¬ 
age  reconstruction,  we  considered  tbis  optical 
property  known  to  estimate  the  best  possible  per¬ 
formance  for  inverse  image  reconstruction  under 
physiological  conditions.  The  two-dimensional  re¬ 
constructed  spatial  map  of  the  fluorescence  yield, 
and  lifetime,  t,  are  shown  in  Figs.  8(a)  and 
8(b),  respectively.  The  image  quality  is  slightly 
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inferior  with  respect  to  both  the  size  as  well  as  the 
shape  of  the  hidden  object.  As  shown  in  Table  3, 
the  mean  value  of  location  of  the  object  according  to 
our  criterion  based  on  T)(j,a>  occurred  at  position 
(59,  58),  consistent  with  the  conditions  used  to  sim¬ 
ulate  the  experimental  data.  The  dimension  of  the 
heterogeneity  based  on  our  arbitrary  definition  for 
identification  (all  grid  points  with  higher 

than  35%  of  the  maximum)  was  683  mm2,  which  is 
close  to  that  used  to  generate  our  simulated  exper¬ 
imental  data.  The  average  values  of  T]|i„  and  t 
in  the  grid  points  that  occupy  the  simulated  object 
converge  within  50  iterations  to  the  values  of  r\jxf,  = 
2.1  X  10“3  mm-1  and  t  =  4  ns,  which  are  slightly 
higher  than  the  values  used  to  generate  the  simulated 
experimental  data  (see  Table  4).  The  average  values 
of  T||jLn  _  and  t  in  the  grid  points  that  occupy  the  sim¬ 
ulated  background  converge  within  50  iterations  to 
values  similar  to  those  reported  in  Case  5. A. 

C.  Two  Heterogeneities  with  Optical  Contrast  from  (t); 
and  (T)|xa^)y  with  No  Absorption  of  Excitation  and 
Fluorescent  Light  by  Chromophores 
In  Case  5.C,  the  same  optical  parameters  were  used 
in  forward  calculations  as  described  in  Subsection 
5. A,  except  that  the  fluorescence  yields,  injjg  for 
objects  1  and  2  were  chosen  as  1  X  10"3  mm-1  and 
2  X  10  3  mm-1,  respectively,  and  lifetimes,  t,  for 
the  objects  were  chosen  as  1  ns  and  2  ns,  respec¬ 
tively.  Two  20-mm-diameter  circular  objects  were 
placed  along  a  diagonal  at  the  coordinates  shown  in 
Table  3  within  a  100-mm-diameter  circular  tissue 
phantom. 

Again  during  the  reconstruction,  no  a  priori  knowl¬ 
edge  of  either  the  object  location  or  the  background 
fluorescence  properties  was  assumed  and  the  values 
of  yield  r\\x„  m  and  lifetime  t  were  found  at  all  the  grid 
points  on  a  33  X  33  grid.  The  two-dimensional  re¬ 
constructed  spatial  map  of  the  fluorescence  yield, 
T||xa  _  ,  is  shown  in  Fig.  9.  The  object  locations  (x,  y) 
obtained  are  given  in  Table  3  and  match  well  with  the 
conditions  used  to  generate  the  simulated  experimen¬ 
tal  data.  The  areas  of  the  objects  from  the  recon¬ 
structed  image  (all  grid  points  with  up, higher 
than  35%  of  the  maximum)  were  381  mm2  (top 
left,  object  1)  and  342  mm2  (bottom  right,  object  2), 
slightly  larger  than  inputs  to  the  forward  problem. 
The  average  values  of  r|p.a  and  t  in  the  grid  points 

that  occupy  the  simulated  object  converge  within  100 
iterations  to  the  values  of  r||ia  =  2  X  10  3  mm-1 
and  =  1.8  X  10"3  mm-1  for  objects  1 

and  2,  respectively,  which  again  is  close  to  the  expected 
values  of  1  X  10"3  mm-1  and  2  X  10"3  mm”1.  Life¬ 
times  of  the  two  objects  were  found  to  be  4. 1  and  3.5  ns, 
respectively,  whereas  the  expected  values  were  1  and  2 
ns,  respectively.  The  quantitative  values  of  t](xq  _ 
and  t  obtained  from  the  inverse  solution  are  currently 
unsatisfactory,  and  research  is  in  progress  to  improve 
the  solution  procedure.  The  background  value  of 
ti|xc,  agrees  well  with  the  expected  value  of  1.0  X 


Fig.  10.  Reconstructed  spatial  map  of  fluorescence  yield,  T|p.aj^m, 
on  a  two-dimensional,  33  X  33  grid  for  the  case  described  in  Sub¬ 
section  5.D.  p.,2i  _  at  the  excitation  wavelength  not  accounting  for 


wavelength  is  the  same.  The  average  value  of  __  within  the 
object  was  2  X  10-4  mm-1  (expected)  and  7.0  X  10~4  mm-1  (re¬ 
constructed).  Spurious  unphysically  high  values  of  have 

been  replaced  by  the  average  background  fluorescence  yield  ob¬ 
tained  from  the  inversion. 


10  5  mm  \  whereas  the  value  of  the  background  life¬ 
time  was  poorly  reconstructed. 

D.  Single  Heterogeneity  with  Optical  Contrast  from  (t), 
and  (T||xa^m)y  with  Absorption  of  Excitation  and  Fluorescent 
Light  by  Chromophores  and  an  Uptake  Ratio  of  20 
The  same  hidden  object  location  as  well  as  optical 
parameters  were  used  as  described  in  Subsection  5. A, 
except  that  a  uniform  background  chromophore  ab¬ 
sorption  coefficient,  |xa^,  of  1  X  10"3  mm-1  as  well  as 
|xa;  at  the  emission  wavelength  of  1  X  10~8  mm"1 
were  used  to  generate  the  simulated  experimental 
data,  and  Gaussian  noise  was  added  as  discussed 
above.  Also,  the  values  of  r||xa  __  for  the  background 
and  the  object  were  chosen  as  1  X  10”5  mm"1  and  2  X 
10"4  mm"1,  thus  giving  an  uptake  ratio  of  20. 
Again,  although  excitation  light  propagation  was  not 
employed  for  image  reconstruction,  we  considered 
this  optical  property  to  be  known  in  order  to  estimate 
the  best  possible  performance  for  inverse  image  re¬ 
construction  under  physiological  conditions.  The 
two-dimensional  reconstructed  spatial  map  of  the  flu¬ 
orescence  yield,  iq|xa  ,  is  shown  in  Fig.  10.  The 
image  quality  with  respect  to  both  the  size  as  well  as 
the  shape  of  the  hidden  object  is  good.  As  shown  in 
Table  3,  the  mean  value  of  location  of  the  object  ac¬ 
cording  to  our  criterion  based  on  T||xa  occurred  at 
position  (61,  57),  consistent  with  the  conditions  used 
to  simulate  the  experimental  data.  The  dimension 
of  the  heterogeneity  based  on  our  arbitrary  defini¬ 
tion  for  identification  (all  grid  points  with  r|pa 
higher  than  35%  of  the  maximum)  was  693  nmt; 
which  is  close  to  that  used  to  generate  our  simu¬ 
lated  experimental  data.  The  average  values  of 
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T)(jLa  and  t  in  the  grid  points  that  occupy  the  sim¬ 
ulated  object  converge  within  50  iterations  to  the 
values  of  tux0  =  7  X  10-4  mm-1  and  t  =  4  ns, 
which  are  sliglitly  higher  than  the  values  used  to 
generate  the  simulated  experimental  data  (see  Table 
4).  The  average  values  of  T)pa  ^  and  t  in  the  grid 
points  that  occupy  the  simulated  background  con¬ 
verge  within  50  iterations  to  values  similar  to  those 
reported  in  Case  5.A.  Nonetheless,  it  is  shown  that 
image  reconstruction  may  be  successful  even  in  the 
case  of  uptake  ratio  as  small  as  20. 

The  performance  of  the  inverse  imaging  algorithm 
may  be  improved  by  including  additional  information 
obtained  from  multifrequency  measurements  and 
from  excitation  wavelength  measurements.  In  ad¬ 
dition,  the  appropriate  weighting  of  grid  point  contri¬ 
butions  based  on  the  signal  strength  may  also 
improve  reconstruction  results. 

6.  Conclusions 

The  fluorescent  yield  and  lifetime  of  an  endogenous 
fluorophore  may  be  sensitive  to  local  environments, 
providing  specificity  for  contrast  of  diseased  over 
normal  tissues  and  optimal  detection  of  disease  by 
using  optical  techniques.32  Our  simulated  experi¬ 
ments  show  that  it  is  possible  to  reconstruct  the 
fluorescent  yield  and  lifetime  of  embedded  fluoro- 
phores  in  tissue-mimicking  scattering  media  from 
frequency-domain  measurements  of  fluorescent 
phase  shift  and  ac  amplitude  or  amplitude  modula¬ 
tion.  As  we  have  shown  in  Subsection  5.C,  the 
reconstruction  of  lifetime  can  be  problematic  by  using 
the  reemitted  fluorescent  signal  at  one  modulation 
frequency.  This  is  in  agreement  with  O’Leaiy  et 
al.,33  who  resorted  to  the  need  for  measurements  in 
the  absence  of  fluorophore  (background  properties)  as 
well  as  in  its  presence  for  lifetime  reconstruction  as 
suggested  from  multipixel  measurements  for  photon 
migration  imaging.34  Currently,  we  are  experimen¬ 
tally  investigating  the  implementation  of  FLI  by  us¬ 
ing  both  excitation  and  fluorescent  wavelengths  as 
well  as  multifrequency  measurements  to  improve  in¬ 
verse  solutions  that  predict  both  fluorescent  yield  and 
lifetime. 

Appendix  A:  Nomenclature 

c,  velocity  of  light 
D(r),  optical  diffusion  coefficient 

f,  modulation  frequency 

I,  identity  matrix 

i,  detector  number,  i  =  1,  20 

j,  Jacobian  matrix  relating  the  sensitivity  of  op¬ 
tical  parameters  to  the  detector  response 

j,  grid  point  number 

jtj,  individual  elements  of  Jacobian  matrix  J 

k,  source  number,  k  =  1,4 

M,  log  of  ac  amplitude  of  modulated  fluorescent 
light 

r,  position 

S(r,  u>),  source  term  for  the  modulated  light  at  position 
r  and  frequency  w 


Greek 

X2,  merit  function  representing  the  least- 
squares  error 

4>(r,  o>),  complex  number  representing  the  photon 
flux  in  the  frequency  domain  at  position  r 
and  frequency  w 

x|,  quantum  yield  of  the  fluorescent  probe  or 
dye 

p,a,  average  absorption  coefficient  of  the  tissue 
p,a  ,  absorption  coefficient  of  the  fluorescence 
light  by  both  the  nonfluorescing  chro- 
mophores  and  fluorophore 
p,a  ,  absorption  coefficient  of  the  excitation  light 
by  both  the  nonfluorescing  chromophores 
and  fluorophore 

|ia  ,  absorption  coefficient  of  the  fluorescence 
light  by  both  the  nonfluorescing  chro¬ 
mophores 

|xa  _  ,  absorption  coefficient  of  the  excitation  light 
by  the  nonfluorescing  chromophores 

^  ,  absorption  coefficient  of  the  excitation  light 
by  the  fluorophore 

|xs. ' ,  effective  scattering  coefficient  of  the  tissue 

0,  phase  shift  of  the  modulated  light  wave 
with  respect  to  the  modulated  wave  at  the 
source 

ct,  standard  deviation  of  the  Gaussian  noise 
representing  the  experimental  uncertainty 
T(r),  lifetime  of  the  activated  probe  or  dye  at 
location  r 

co,  angular  modulation  frequency,  given  by 
2tt f 

Subscripts 

obs,  observed  or  experimental  data 
x,  excitation  light 
m,  fluorescence  or  emission  light 

Appendix  B:  Jacobian  Matrices 

We  provide  more  details  about  the  elements  of  the 
Jacobian  matrices  introduced  in  Section  4.  Element 
jij  of  a  Jacobian  matrix  represents  the  sensitivity  of 
the  response  of  detector  i  to  changes  in  the  optical 
property  at  grid  point  j.  Here  as  an  example  we  show 
the  response  of  detector  16  for  source  A  (see  Fig.  1)  to 
changes  in  t  and  Tip,  at  all  the  grid  points  for  the 
case  described  in  Subsection  5.A  for  the  first  iteration 
during  inversion  when  |xt;  was  chosen  as  zero.  Sim¬ 
ilar  results  are  observed  for  other  cases.  Figures 
11(a)  and  ll(b)=show  the  elements  of  Jacobian  matri¬ 
ces  J(0,  t)  and  J(M,  t),  respectively,  where  t  at  each  of 
the  grid  points  was  increased  by  5%.  Similarly,  =Fig. 
11(c)  shows  the  elements  of  Jacobian  matrix  J(M, 
•qp.^  ),  where  (rnx^^J  at  each  of  the  grid  points  was 
increased  by  1%.  Most  of  the  elements  of  J(0,  t) 
and  J(Af,=T|p,a^J  are  positive,  whereas  most  of  the  ele¬ 
ments  of  J (M,  t)  are  negative.  This  is  according  to  our 
expectation  of  systems  with  no  scattering.  The  de¬ 
pendence  on  scatter  has,  of  course,  been  taken  into 
account  in  the  above  example. 
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Fig.  11.  Jacobian  (a)  j(9,  t),  (b)  J(M,  t),  (c)  J(M,  for  the 

case  described  in  Subsection  5.A  for  source  A,  detector  16,  and 
iteration  1.  During  the  computation  of  the  Jacobians,  the  values 
of  t  for  (a)  and  (b)  and  tip,a  for  (c)  at  each  individual  grid  point 
were  increased  by  5%,  5%,  and  1%,  respectively. 
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Abstract 

From  analytical  and  numerical  solutions  predicting  the  scattering  of 
diffuse  photon  density  waves  and  from  experimental  measurements  of  changes  in 
phase-shift,  6,  and  AC  amplitude  demodulation,  M,  due  to  the  presence  of  single 
and  double  cylindrical  heterogeneities,  we  show  that  second  and  higher  order 
perturbations  can  impact  the  prediction  of  the  propagation  characteristics  of 
diffuse  photon  density  waves.  Our  experimental  results  for  perfect  absorbers  in  a 
lossless  medium  suggest  that  the  performance  of  fast  inverse  imaging  algorithms 
which  employ  first  order  Born  or  Rytov  approximations  may  have  inherent 
limitations  as  compared  to  inverse  solutions  that  employ  iterative  solutions  of  a 
linear  perturbation  equation  or  numerical  solutions  of  the  diffusion  equation. 

*To  whom  correspondence  should  be  addressed 
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1.0  Introduction 


With  the  development  of  near-infrared  red  (NIR)  emitting  laser  diodes 
and  the  realization  that  NIR  light  can  travel  several  centimeters  through 
tissue,  numerous  groups  have  embarked  upon  the  development  of  optical 
imaging  as  a  new  medical  imaging  modality.  While  approaches  vary  from 
monitoring  the  vanishingly  small  component  of  coherent  light  using 
sophisticated  techniques  of  Kerr-filtering,1  time-gating,2  and  conservation  of 
light  polarization,3  other  approaches  focus  on  monitoring  the  predominant 
optical  signal  re-emitted  from  tissues:  the  multiply  scattered  signal. 
Continuous  wave  (CW),4  time-domain,5  and  frequency-domain6'8 
measurements  of  multiply  scattered  light  have  been  performed  in  simulated 
and  experimental  tissue  phantom  studies  as  well  as  in  human  studies.9 
However,  the  use  of  these  measurements  to  reconstruct  maps  of  internal 
tissue  optical  properties  for  diagnostic  imaging  has  been  problematic,  since  the 
geometric  correlation  between  incident  and  detected  radiation  is  destroyed  by 
multiple  scattering. 

While  some  investigators  have  used  direct  image  reconstruction 
approaches  employing  measured  optical  data  to  directly  form  an  image,9'11 
others  have  sought  full  solution  to  the  inverse  imaging  problem,  which  relates 
external  time-dependent  measurements  made  at  the  periphery  to  an  optical 
property  map  of  interior  volumes  via  the  optical  diffusion  equation6'8* 12'16. 
Two  approaches  to  solve  the  inverse  solution  have  been  adopted.  In  the  first 
approach,  the  fluence  associated  with  a  propagating  photon  density  wave 
launched  at  source  position,  ps  and  detected  at  the  tissue  periphery  at  position 
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Pd  (denoted  <3>(ps  ,pd)),  is  related  to  the  fluence  assumed  in  the  absence  of  any 
optical  heterogeneities  (denoted  OinC(ps  ,pd))  and  the  internal  optical  property 
perturbation  map,  A(p),  through  an  Fredholm  integral  equation  of  the  first 
kind: 

<D(ps,pd)  =  Oinc(ps,pd)  + jG(p,pd)O0(ps,pd)A(p)a3p  (1) 

where  G(p,pd)  is  the  Green's  function  to  the  diffusion  equation 

(=exp  ^  pa  +  ito/cn)/D^_  describing  the  propagation  of  light  from 
[  47uD(p  -pd)  j 

position  p  to  the  detector  at  pd  and  <J>0  is  the  incident  wave  impinging  at 
position  p.  If  one  assumes  that  the  incident  wave  impinging  upon  the  position 
p  can  be  approximated  by  the  fluence  predicted  in  the  absence  of 
heterogeneities  (i.e.,  <3>0  — >  known  as  the  Born  approximation),  then  upon 

measurement  of  fluence  at  a  variety  of  source-detector  separations  and  upon 
discretization  of  Eqn  (1),  an  optical  map  of  perturbations,  A(p),  can  be 

obtained  to  a  first  order  approximation.15  However,  the  Born  approximation 
<I>0  — »  Oinc  does  not  account  for  second  and  higher  order  effects  which  may 

arise  from  the  re-scattering  of  photon  density  waves  associated  with 
neighboring  inhomogeneities.  In  addition,  this  approach  assumes  that 
perturbations  in  optical  properties  at  a  position  pj  does  not  influence  the 
propagation  of  light  from  position  pj+i  to  detector  pd  as  described  by  the 
Green's  function,  G(p,pd). 

The  Born  iterative  method  (BIM)  tends  to  account  for  strong 
perturbations  and  for  second  and  higher  order  scattering  effects  by  using  the 
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first  order  approximation  of  optical  property  perturbations,  A(p),  to  recalculate 
Oo  or  the  fluence  incident  at  position  p.  Yao  and  coworkers  have  shown  that 
the  Born  iterative  method  tends  to  compensate  for  the  underprediction  of  the 
scattering  and  absorption  properties  of  a  single  heterogeneity  in  an  otherwise 
uniform  medium  that  would  occur  when  the  Born  approximation  (or  single 
iteration)  is  used.  In  addition  to  the  Bom  iterative  methods,  the  Distorted  Bom 
iterative  (DIBM)  method  represents  a  refinement  in  that  it  also  recompiles  the 
Green's  functions,  G(p,pd)  to  reflect  changes  in  light  propagation  and  should 
speed  convergence.  Yet  to  date,  there  has  been  no  investigation  to  address  how 
these  iterative  reconstruction  algorithms  impact  the  resolution,  or  more 
concisely  how  these  linear  perturbation  approaches  can  be  used  to  accurately 
image  closely  positioned,  multiple  heterogeneities  between  which  non-linear, 
second  and  higher  order  scattering  effects  exist. 

In  contrast  to  this  inversion  approach,  Jiang,  et  aZ.7>8  and  Arridge,  et 
al., 14  have  utilized  full  numerical  solution  to  the  diffusion  equation  which 
describes  the  interdependence  of  voxel  optical  properties  and  their  contribution 
to  the  re-emitted  optical  signal  detected  from  CW  and  frequency-domain 
measurements.  Using  this  approach,  the  second  and  higher  order  scattering 
arising  from  multiple  inhomogeneities  which  are  not  accounted  for  in  the  first 
iterations  of  Bom  and  Rytov  approximations,  are  incorporated  in  the  forward 
and  inverse  solutions.  The  inversion  consists  of  relating  spatial  changes  of 
optical  properties  on  detected  frequency-domain  measurements  through 
numerical  solution  of  the  diffusion  equation  and  then  solving  an  update  to  the 
map  of  optical  properties  from  the  difference  between  the  signals  that  are 
measured  and  those  that  are  predicted  by  the  forward  solution  to  the  diffusion 
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equation.  It  is  important  to  note  that  convergence  upon  the  optical  properties 
is  achieved  with  these  numerical  approaches  and  not  with  the  iterative  Born 
and  Rytov  approaches.  Nonetheless,  while  the  computational  investment  of 
iterative,  but  analytical-based  inversions  is  less  than  that  involving  full 
numerical  solution  of  the  diffusion  equation,  the  relative  performance  of  these 
two  inverse  solution  approaches  has  yet  to  be  evaluated. 

For  this  reason,  we  embarked  on  a  study  to  assess  the  contributions  of 
neighboring  heterogeneities  using  experimental,  numerical,  and  analytical 
computations  of  scattered  photon  density  waves  from  perfect  light  absorbing 
cylinders.  Specifically,  we  experimentally  and  computationally  monitor  the 
multiple  scattering  of  photon  density  waves  between  two  neighboring  perfect 
cylindrical  absorbers  embedded  in  a  tissue-mimicking  scattering  medium  in 
order  to  assess  the  higher  order  perturbation  effects  upon  a  re-emitted, 
detected  photon  density  wave.  In  the  following,  we  briefly  review  the  theory  of 
higher  order  perturbation  effects  as  predicted  from  the  experimental 
frequency-domain  measurements  as  well  as  from  the  Helmholtz  equations.  In 
addition,  we  present  experimental  measurements  and  numerical  solutions  of 
the  optical  diffusion  equation  which  show  the  errors  introduced  by  neglecting 
second  order  effects  which  can  be  significant.  These  errors  can  define  the 
limits  of  resolution  for  two  perfect  absorbers  in  the  inverse  solution  approaches 
which  do  not  employ  BIM,  DBIM,  or  the  full  numerical  solution  of  the  optical 
diffusion  equation.  We  comment  upon  these  effects  when  contrast  is  owing  to 
mechanisms  other  than  a  perfect  absorber. 
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2.0  Theoretical  background  for  perturbations  associated  with  diffuse 
photon  density  waves 

Our  work  to  assess  the  contributions  of  higher  order  perturbations  has 
been  motivated  by  the  work  of  Boas,  et  aZ.15  Their  approach  for  image 
reconstruction  from  diffusely  propagating  photon  density  waves  employs 
analytical  expressions  to  describe  the  complex  fluence  of  a  propagating  photon 
density  wave,  <5>(p),  in  the  presence  of  m  heterogeneities  by  superposition  of 
incident,  <&inc(p),  and  scattered  waves,  <E>gCat  k  (p)  ,  from  body  k  and  of  order  n: 

oo  m 

<D(p)  =  (p)  +  I  I  <l£cat  k  (p)  (2) 

n=lk=l 

Figure  1  is  reproduced  from  the  work  of  Boas  and  coworkers  to  pictorially 
describe  the  second  order  effects  due  to  the  scattering  of  a  photon  density  wave 
between  two  objects  and  to  delineate  their  contribution  to  the  detected  photon 
density  wave,  <J>(p).  Second  and  higher  order  scattering  effects  arise  from  the 
re-scattering  of  an  incident  wave  between  the  two  bodies.  These  authors 
assume  that  second  order  scattering  effects  (denoted  by  the  dotted  line  in 
Figure  1)  are  negligible  under  most  circumstances.  In  our  work,  we  explore  this 
assumption  using  experimental,  numerical,  and  analytical  predictions  of 
second  order  contributions. 

2.1  Analysis  of  experimental  measurements  of  photon  density  waves 

In  order  to  measure  second  order  contributions,  we  experimentally 
measured  the  detected  photon  density  wave  in  the  presence  of  one  light 
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absorbing  object  alone,  &kl>  object  two  alone,  0>k2\  and  with  both  objects 
present,  <5>kl,k2  as  a  function  of  object  positions,  pi;  and  p2,  and  separation 
between  the  two  objects,  pi-p2  (See  Figure  2).  From  Eqn  (2),  expressions  can 
be  written  for  measurements  of  <£>kl  and  d>&2  made  at  detector  position,  pd,  in 
the  presence  of  single  objects: 


$ki(pd)  =  ^inc(Pd)  +  ^scat,kl(pd)  (3) 

$k2  (Pd )  =  ®  Inc  (Pd  >  +  <cat,k2  <Pd  >  (4) 

where  kl  denotes  the  presence  of  the  first  body  alone  and  k2  denotes  the 
presence  of  the  second  body.  With  simple  algebraic  manipulation  of  Eqns.  (2) 
through  (4),  an  expression  accounting  for  higher  order  scattering  effects  can  be 
written  for  the  photon  density  wave  in  the  presence  of  both  objects, 
®*l,*2(pd): 


®kl,k2  (Pd )  =  ^kl(Pd  )  +  ®k2  (Pd )  “  °inc  (Pd  )  + 

(5) 

^>scat,kl,k2^Pd^ +  ^scat.k  l)k2^Pd  ^  ‘ ' 

For  the  purposes  of  this  study,  it  is  assumed  that  third  and  higher  order 
scattered  waves  are  considered  to  have  insignificant  contributions  to  the 
measured  fluence  when  compared  to  first  and  second  order  scattered  waves. 

From  frequency-domain  measurements  of  phase-shift,  0,  and  AC 
amplitude,  M,  values  of  the  complex  fluence,  O=Me10  can  be  obtained  (i)  in  the 
presence  of  the  first  object  alone  Q>kl( Pd)>  (ii)  in  presence  of  the  second  object 
alone  <J>&2(Pd)>  (iii)  in  the  presence  of  both  objects,  ®kl,k2( Pd)>  and  in  the 
absence  of  any  inhomogeneities,  OincCpd)-  If  second  order  effects,  i.e., 
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^scat  ki,k2^Pd)’  are  negligible,  then  from  Eqn  (5)  it  follows  that  measurements 
of  ®kl,k2( Pd)  should  be  predicted  by  measurements  of  <£jfei(pd),  ^^(pdX  and 
^incCPd): 


$kl,k2  (Pd )  =  ®k=l(Pd )  +  $k=2  (Pd )  -  ®inc  (Pd  )  (6) 


Since  we  report  our  results  in  terms  of  phase-shift  and  amplitude 
demodulation,  Eqn  (6)  can  be  written  as: 


ekl,k2(Pd)  =  tan  1 


Mkisin9kl  +  Mk2  sinek2  -  sin  9^ 
Mki  cos  9k  i  +  Mk2  cos  9k2  -  cos  9inc  _ 


(7) 


Mkl,k2(Pd)  = 


1 


[Mkl  cos  9kl  +  Mk2  cos  9k2  -  cos  9^  r  + 

[Mklsin9kl  +  Mk2  sin9k2  -  Minc  sin9inc]J 


(8) 


2.2  Analysis  of  analytical  predictions  of  phase-shift  and  amplitude 
modulation  owing  to  two  perfect  absorbing  cylinders 

The  complex  fluence  describing  the  propagation  of  a  photon  density 
wave  can  also  be  computed  analytically  from  the  Helmholtz  equations.15*16 
The  complex  fluence  propagating  from  a  source  point  at  ps,  through  an  infinite 
random  medium,  and  received  at  position,  p,  in  the  absence  of  an  optical 
heterogeneity  is  given  by: 
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^inc^P)  —  ^source (Ps) 


exptt-j  ^at^Z^n[p  _pd]} 
4tiD[P  -pd] 


where  Ss0urce  describes  the  phase  and  strength  of  modulation  of  the  source 
located  at  position  ps;  cn  is  the  speed  of  light  in  the  medium;  and  D  is  the 
"optical  diffusion  coefficient"  which  is  governed  by  the  absorption,  pa,  and 
isotropic  scattering,  p's,  coefficients  of  the  continuous  or  homogeneous  medium: 

D  =  -7-  1  ~n  do) 

3[pa+Ps] 

In  this  study,  we  approximate  a  point  modulated  source  at  the  periphery 
of  a  large  cylinder  as  a  point  modulated  source  located  at  the  interface  of  a 
semi-infinite  medium.  In  this  case,  the  fluence  is  written  in  cylindrical 
coordinates  and  at  the  surface  (z=0)  given  by:17 


1 

M-a+ico/cn  / 

1 

D  \ 

-((Ps-Pd)2+(z-z0)2)2 


^inc  (Pd )  ~  ^source  (Ps ) ' 


/((Ps-Pd)2+(Z-Z0)2 


Pa+ico/cn  / 

1 

D  \ 

((ps  -  Pd  )2  +  (z  +  2o )2  )2 


[  #s-Pd)2+(z  +  z0)2 

where  z0  is  one  isotropic  scattering  length  (l/|x's). 


The  first  order  (n=l)  scattered  wave  from  the  kl  infinitely  long  cylinder 
in  an  infinite  medium  is  given  by:15>18 
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°scatkl(Pd)  =  -°inc  I  Jcos(mfi)cos(pz)K 

m=l  0 


m 


2  .  r— Ha  +  fo/Cnl  1 


P  + 


D 


Pd 


>K 


m 


p2+r_Ha±K»Z£n  |Ps 


DxIm(x)Im(y)-D|,1yIm(x)Im(y) 
LDxKm(x)Im(y)-DklyKm(x)Im(y)J 


dp 


(12) 


where  is  the  optical  diffusion  coefficient  within  the  cylinder,  k;  Im  and  Km 
are  modified  Bessel  functions;  x  =  -p2+  -  P,a_+*(D^cn  a^. 


y  = 


n2  + 

P  + 

L  Dk  J 

a^;  and  ak  is  the  radius  of  cylinder  kl.  Radius  pd, 


angle  $ ,  and  length  z  denote  the  coordinates  of  the  detector  relative  to  the 
center  of  the  infinite  cylinder,  kl  (Figure  2).  The  incident  wave  upon  the  infinite 
cylinder,  <X>inc ,  is  computed  from  Eqn  (11).  A  similar  expression  can  be  written 
for  0?cat,k2' 


Second  order  scattering  contributions  were  calculated  by  evaluating  the 
scattered  wave  originating  from  the  second  object  ( k2 )  as  the  incident  wave 
upon  the  first  object  (kl).  In  other  words, 
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0scat,kl(Pd)  =  -^scat,k2  5  Jcos(mfl)cos(pz)Kmf Jp2  +  f-^^/Cnlpd 

m=l  0  U  v  u  )  ) 


.  D  x  Tm  (x)  lm  (y)  ~  Dkl  y  ^  (x)  4  (y)  d 
#|_DxKm(X)Im(y)-DklyKm(x)Im(y)J  P 


where  the  incident  wave  upon  cylinder  kl  is  now  the  first  order  scattered  wave 
<&scat  k2  that  is  computed  from  Eqn  (12).  A  similar  expression  can  be  written 

for  *^scat,k2‘  Since  the  optical  properties  and  diameter  of  the  two  cylinders 

were  identical  in  this  study,  we  consider  the  case  where  =  D2  and  ai=a2- 


Using  the  approach  described  above,  the  fluence  detected  at  position  pd 

can  be  computed  inclusive  of  first  order  scattering  effects  (i.e., 

2 

$(pd)  =  ^inc  (Pd ) +  X  ^scatk^Pd^*  as  weU  as  second  order  scattering  effects 
k=l 

22 

(i.e.,  <J>(pd)  =  d>inc(pd)+  X  X  ^scatk^d))-  From  final  values  of  complex 

n=l  k=l 

fluence,  the  phase-shift  and  amplitude  demodulation  can  be  predicted  from  the 
simple  relationships: 


0(pd)  =  tan 


-1  Im[P(pd)] 
Re[<t>(pd)j 


M(Pd)  =  -y[Im{<l)(pd)}]2  +[Re{0(pd)}]5 
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3.0  Materials  and  Methods 


3.1  Experimental  measurements  of  phase-shift  owing  to  two 
cylindrical  absorbers 

In  order  to  experimentally  measure  M&j^Pd),  M^pd),  M&2(pd)> 
Minc(pd),  ®kl,k2( Pd),  0Ai2(pdX  0fc2(pd),  and  0inc(pd)>  frequency  domain 
measurements  were  made  using  an  apparatus  employing  picosecond  pulsed 
light  at  780  nm  with  an  average  power  of  1.3  watts.  Details  of  the  apparatus 
are  described  elsewhere.19  The  phantom  consisted  of  a  plexiglass  cylinder 
(with  a  16.5  cm  diameter  and  20  cm  height),  filled  with  a  0.5%  Intralipid 
solution  (Kabi  Pharmacia,  Inc.,  Clayton,  NC).  As  illustrated  in  Figure  3,  light 
was  delivered  to  a  peripheral  point  on  the  cylinder  using  a  1000  micron  fiber 
(HCP-M1000T-08,  Spectron  Specialty  Optics,  Co.,  Avon,  Conn)  and  collected 
via  a  second  fiber  located  4  circumferential  cm  from  the  incident  source. 
Bakelite  plastic  rods  (diameter  =  3.175  mm)  were  painted  black  to  provide 
perfectly  light  absorbing  inhomogeneities.  Measurements  of  0(pd)  and  M(pd) 
at  80  MHz  were  conducted  as  the  rods  were  moved  in  tandem  along  the  plane 
perpendicular  to  a  line  connecting  the  source  and  detector.  Positioning  was 
achieved  with  a  motion  controller  (PMC200-P,  Newport  Corp.,  Irvine,  CA)  and 
a  motorized  actuator  (Newport  850B).  The  actuator  position  was  accurate  to 
within  0.0050  mm.  Phase  and  AC  modulation  were  recorded  as  the  distance 
between  the  center  of  the  first  absorbing  cylinder  and  the  wall  of  the  phantom 
was  varied  from  0  cm  to  5  cm  in  20  increments.  When  the  two  objects  were 
moved  in  tandem,  their  separation  distances  (pi  -  P2)  between  the  centers  of 
the  two  perfect  absorbers  were  6  mm,  10  mm,  and  20  mm.  Three 
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measurements  were  taken  at  each  position.  Measurements  of  phase-shift  and 
AC  amplitude  demodulation  are  reported  relative  to  the  absence  case  or  0inC  or 
Mjnc- 

3.2  Analytical  prediction  of  phase-shift  and  amplitude  modulation 
owing  to  two  absorbing  cylinders 

In  addition  to  the  measurements  of  second-order  interactions  with  the 
experimental  approach  described  above,  predictions  of  interactions  between 
absorbers  were  computed  analytically  using  Eqns.  (11)  through  (15)  upon 
employing  a  modified  version  of  an  algorithm  written  by  D.A.  Boas  and 
available  through  the  internet  (http://www.lrsm.upenn.edu/pmi/PMI/pmi/html). 
Complex  fluence  incorporating  first  and  second  order  scattering  contributions 
were  used  to  calculate  phase-shift  owing  to  two  infinite  cylinders  that 
effectively  act  as  perfect  absorbers  (|ia  =  2  cm'1,  (is  =10  cm  )  in  a  turbid, 
semi-infinite  medium  mimicking  our  phantom  (|ia  =  0.02  cm'1,  p.s'=10  cm,-1). 
Zero  fluence  boundary  conditions  were  assumed  and  employed  in  the  algorithm. 

3.3  Numerical  prediction  of  phase-shift  and  amplitude  demodulation 
owing  to  two  absorbing  heterogeneities 

In  addition  to  experimental  measurements  and  analytical  calculations, 
two  dimensional  finite  element  computations  were  conducted  in  order  to  predict 
the  phase-shift  and  amplitude  demodulation  owing  to  first  order  and  second 
order  interactions.  The  computations  were  performed  on  an  Ultra  2  Sun  Sparc 
workstation  using  a  MATLAB®  partial  differential  equation  toolbox.  The 

Sevick-Muraca 
page  13 
May  20,  1997 


frequency  domain  diffusion  equation  for  light  propagation  was  solved  to 
determine  the  fluence  from  an  80  MHz  sinusoidally  modulated  light  source. 
The  simulated  phantom  was  an  infinite  cylinder  16.5  cm  in  diameter.  The 
optical  properties  of  the  medium  were  set  to  mimic  the  experimental  conditions 
of  0.02  cm*l  for  the  absorption  coefficient,  |ia>  and  10  cm'l  for  the  isotropic 
scattering  coefficient,  |i's.  The  nearly  perfect  absorbers  in  these  computations 
were  modeled  as  infinite  cylindrical  rods  and  approximated  by  increasing  the 
absorption  coefficient  to  one  hundred  times  that  of  the  surrounding  medium  (|ia 
=  2.0  cm'l).  The  scattering  coefficient  for  the  object  was  set  to  that  of  the 
surroundings  (p's  =  10  cm'l)  to  mimic  the  analytical  computations.  The 
phantom  was  discretized  into  66048  triangular  elements  containing  33281 
nodes.  A  partial  current  boundary  condition  was  used  to  approximate  light 
reflection  and  transmission  across  the  boundary  of  the  phantom  as  would  be 
expected  in  the  experimental  measurements.  While  our  meshing  did  not  permit 
representation  of  the  perfect  absorber  as  a  volume  excluded  for  light  transport, 
we  approximated  the  heterogeneity  with  a  high  absorption  coefficient  since 
others  have  shown  that  the  propagation  characteristics  are  comparable.18  It 
is  noteworthy  that  while  these  computations  do  not  exactly  mimic  the 
experimental  measurements  described  below,  they  nonetheless  adequately 
describe  the  contributions  of  second  order  scattering  effects.  The  forward 
solution  was  obtained  for  each  absorbing  cylinder  alone  and  then  in 
combination.  The  results  were  analyzed  similar  to  the  experimental  results 
using  Eqns.  (7)  and  (8). 


Sevick-Muraca 
page  14 
May  20,  1997 


4.0  Results  and  Discussion 

4.1  Experimental  results 

Figures  4a  through  4c  and  Figures  5a  through  5c  illustrate  the  phase- 
shift  and  amplitude  demodulation  changes  as  a  function  of  position  of  the  pair 
of  perfect  absorbers  whose  centers  are  separated  by  6  mm,  10  mm,  and  20 
mm.  The  symbols  denote  individual  measurements  of  phase-shift  difference  in 
the  presence  of  the  two  hidden  objects  while  the  symbols  connected  by  the  lines 
denote  the  phase-shift  and  amplitude  demodulation  calculated  from  Eqns  (7) 
and  (8)  with  measurements  made  in  the  absence  and  in  the  presence  of  a  single 
individual  absorber.  The  measured  values  of  phase-shift  and  amplitude 
demodulation  denoted  by  the  open  symbols  are  therefore  reflective  of  higher 
order  contributions  while  the  solid  line  reflects  only  first  order  contributions. 
Paired  student  t-test  shows  that  there  is  significant  difference  (p<0.005) 
between  the  set  of  phase-shift  and  amplitude  demodulation  measurements  and 
that  predicted  by  Eqns  (7)  and  (8)  for  two  perfect  light  absorbing  objects 
separated  by  6, 10,  and  20  mm  and  positioned  at  varying  distances  away  from 
the  source  and  detector  (Figures  4a  through  5c).  Our  results  also  show  that 
second  order  effects  are  greatest  for  two  absorbing  heterogeneities  spaced  6 
mm  apart  and  become  smaller  in  magnitude  at  10  and  20  mm.  Paired  t-test 
indicate  there  are  significant  differences  at  the  99.5%  confidence  levels 
(indicated  by  the  astericks  in  Figure  4)  between  individual  experimental  phase- 
shift  values  that  reflect  first  and  higher  order  scattering  contributions  to  the 
detected  signal  and  those  computed  values  that  are  indicative  of  first  order 
effects  only.  From  the  data  for  6  mm  absorber  separation  illustrated  in  Figure 
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4a,  it  can  be  seen  that  the  actual  experimental  measurements  of  phase-shift 
change  due  to  two  objects  is  smaller  than  that  predicted  by  Eqn  (7)  in  which 
second  and  higher  order  perturbations  are  not  accounted  for.  At  greater 
distances  away  from  the  source  and  detector,  agreement  between 
experimental  and  predicted  phase-shift  change  suggests  that  second  order 
effects  may  indeed  be  negligible,  even  in  the  case  of  the  smallest  separation, 
but  only  in  a  region  where  the  objects'  contribution  to  the  detected  signal  is 
comparatively  small.  This  is  reasonable  since  second  order  effects  are 
expected  to  increase  with  proximity  to  the  source  and  detector.  From  our 
studies,  it  appears  that  separations  greater  than  20  mm  are  necessary  for 
their  accurate  resolution  via  perturbative  reconstruction  approaches  when 
perfect  absorbers  are  involved. 

4.2  Analytical  computations 

Our  experimental  results  are  also  validated  by  the  analytical  predictions 
which  account  for  first  and  second  order  scattering  contributions.  Figures  6a 
through  6c  depict  the  qualitative  trends  seen  in  the  experimental  data 
presented  for  phase  difference  in  Figures  4a  through  4c.  There  appears  to  be 
good  agreement  between  the  trends  observed  experimentally  with  varying 
separation  distances  and  analytical  predictions.  Specifically,  higher  order 
contributions  significantly  perturb  the  detected  phase-shift  values  when  the 
separations  between  the  center  of  two  absorbing  cylinders  are  6  mm  (Figure 
6a)  and  10  mm  (Figure  6b).  The  contribution  of  higher  order  scattering  from 
cylindrical  absorbers  is  not  significant  when  the  separation  distance  is  20  mm 
(Figure  6c).  However  upon  comparison  to  experimental  results,  there  are 
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differences  in  the  absolute  magnitude  and  shape  of  the  phase  shift  change 
versus  object  distance.  These  distances  are  most  likely  due  to  the  differences 
in  the  boundary  conditions,  geometry,  and  absorber  strength  between  the 
analytical  and  experimental  results.  Nonetheless,  the  trends  confirm 
experimental  results  that  the  presence  of  second  order  perturbations  reduces 
the  phase  shift  change  owing  to  two  light-absorbing  bodies.  Neglecting  second 
order  effects  could  cause  an  underestimation  of  absorption  strength  or  size 
when  reconstruction  algorithms  based  upon  first  order  perturbations  are  used. 

In  addition,  we  investigated  the  variation  of  second  order  perturbations 
with  modulation  frequency  as  shown  in  Figure  7.  The  simulated  data  of  phase- 
shift  change  versus  position  of  the  heterogeneities  is  depicted  for  two 
cylindrical  absorbers  (d  =  3.125  mm,  |ia=2  cm’1 )  separated  by  6  mm  in  a  semi¬ 
infinite  medium  Again,  the  phase-shift  is  reported  relative  to  the  absence  case 
and  the  distance  is  reported  from  the  wall  to  the  center  of  the  first  cylinder. 
Three  frequencies  were  evaluated:  80, 160,  and  240  MHz.  The  absolute  value 
of  phase  shift  increases  with  frequency,  but  the  contribution  of  second  order 
effect  remains  roughly  the  same  absolute  value.  Consequently,  the  error  in 
assuming  negligible  second  order  scattering  effects  becomes  smaller  at 
increasing  frequencies.  This  is  expected  since  increased  damping  of  a  re¬ 
scattered  photon  density  wave  occurs  at  higher  modulation  frequencies.  While 
these  results  intimately  depend  upon  the  choice  of  optical  properties,  they 
nonetheless  point  out  that  the  error  in  neglecting  second  order  scattering 
contributions  in  analytically-based  reconstructions  becomes  smaller  with 
increasing  modulation  frequencies.  Of  course,  the  error  reduction  occurs  at  the 
expense  of  interrogating  a  smaller  volume  of  tissue  with  a  source  modulated  at 
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an  increased  frequency.20 


4.3  Finite  element  computations 

Our  experimental  results  were  not  only  validated  by  analytical 
predictions  but  also  with  finite  element  computations.  Figures  8a  -  9c  show 
the  numerical  computations  that  correspond  to  the  experimental  phase-shift 
data  presented  in  Figures  4a-5c.  Again  there  appears  to  be  good  agreement 
between  the  trends  of  the  experimental  and  analytical  results  with  that  of  the 
numerical  solutions.  Also,  higher  order  contributions  perturb  the  detected 
phase-shift  values  when  the  separation  distance  between  the  two  rods  are  6 
mm  (Figure  8a)  and  10  mm  (Figure  8b).  The  contribution  of  higher  order 
scattering  from  the  cylindrical  absorbers  is  not  significant  when  the  separation 
distance  is  20  mm  (Figure  8c).  Figures  9a-  9c  show  the  amplitude 
demodulation  relative  to  an  absence  condition.  The  modulation  data  also 
shows  similar  results  to  that  obtained  experimentally.  These  computations 
results  confirm  that  the  presence  of  second  order  perturbations  are  important 
for  two  light  absorbing  objects  that  are  less  than  20  mm  apart. 

5.0  Conclusions 

In  summary,  our  experimental  measurements  show  that  the 
contributions  of  higher  order  scattering  of  propagating  photon  density  waves 
may  not  always  be  insignificant.  While  our  analytical  and  numerical 
computations  do  not  exactly  reproduce  experimental  conditions  (i.e.,  2-D  finite 
element,  semi-infinite  geometry,  etc.),  they  nonetheless  demonstrate  that  the 
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experimental  trends  can  be  attributed  to  second  order  effects.  Under 
conditions  of  high  contrast  owing  to  absorption,  analytical  approaches  to  the 
inverse  imaging  algorithm  may  restrict  the  resolution  and  sensitivity  of 
biomedical  optical  imaging  performed  in  the  'frequency-domain.  An  analogy  can 
also  be  drawn  for  time-domain  and  CW  reconstruction  approaches  which  do 
not  deploy  full  solutions  to  the  diffusion  equation  to  account  for  the 
interdependence  of  voxel  optical  properties  on  measured  fluences.  Certainly, 
our  results  are  based  upon  the  worst  case  scenario  of  perfect  absorbers,  and 
may  have  less  impact  on  the  image  reconstructions  involving  imperfect 
absorbers.  Under  conditions  in  which  multiple  heterogeneities  are  contrasted 
from  their  surroundings  on  the  basis  of  scattering  or  fluorescence,  the  first 
order  perturbation  assumption  may  not  be  as  restrictive.  Indeed,  if  the  non¬ 
linearity  associated  with  second  order  scattering  effects  are  small,  then  non¬ 
iterative  Bom  and  Rytov  approximations  are  especially  attractive  since  image 
reconstruction  are  not  computational  intensive.  Nonetheless,  our  results 
suggest  that  inverse  imaging  algorithms  that  depend  solely  upon  first  order 
perturbations  due  to  local  changes  in  tissue  absorption  properties  may  not 
provide  as  accurate  a  reconstruction  when  compared  to  those  algorithms  that 
depend  upon  the  full  numerical  calculation  of  the  diffusion  equation  with 
specified  boundary  conditions. 
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6.0  Nomenclature 


ak 

cn 

D 

Dk’ 

In 

In1 

Kn 

Kn’ 

M 

m 

P 

Ssource 


x 


y 

z 


Radius  of  cylinder  k,  [cm]. 

Speed  of  light  through  the  medium,  [cm/sec]. 

Optical  diffusion  coefficient  of  homogeneous  medium, 
[cm]. 

Diffusion  coefficient  inside  the  cylinder,  [cm]. 
Modified  Bessel  function. 

Derivative  of  modified  Bessel  function,  In- 

Modified  Bessel  function. 

Derivative  of  modified  Bessel  function,  Kn. 

AC  demodulation  of  the  incident  light, 

[mW/cm2] . 

Number  of  objects. 

Integration  variable  in  Helmholtz  equation, 
describing  scatter  from  an  infinite  cylinder,  [cnr1]. 
Strength  of  modulated  source  at  position  ps? 
represented  as  a  complex  number  of  amplitude  and 
phase [mW]  . 

Variable  in  Helmholtz  equation 

describing  scatter  from  an  infinite  cylinder,  [cm'1]. 

variable  in  Helmholtz  equation 

describing  scatter  from  an  infinite  cylinder,  [cm-1]. 

Axial  direction  of  the  cylindrical  heterogeneities, 

[cm]. 
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Greek: 


<E>  Photon  fluence,  [mW/cm2]. 

Ojnc  Photon  fluence  of  incident  wave  or  in  the  absence  of 

heterogeneities,  [mW/cm2]. 

Os  cat  Photon  fluence  arising  from  scattered  wave, 

[mW/cm2]. 

p  Position  vector,  [cm]. 

ps  Position  of  the  source,  [cm], 

pd  Position  of  detector,  [cm], 

pk  Center  position  of  object  k,  [cm], 

0  Phase  shift  of  light  wave,  [degrees  or  radians] . 

■0  Angle  between  ps  and  pd- 

|ia  Absorption  coefficient  of  homogeneous  surroundings, 

[cm-1]. 

Ps'  Isotropic  scattering  of  homogeneous 

surroundings,  [cm-!]. 

Pa'  Absorption  coefficient  of  cylinder,  [cm-!]. 

Ps"  Isotropic  scattering  of  cylinder,  [cm-!]. 

subscripts  and  superscripts 

n  Order  of  perturbation  or  nth  order  scattering  effect, 

k  Index  denoting  cylindrical  object  one,  kl,  two,  k2, 

or  both  objects  kl,k2. 

inc  Index  denoting  absence  measurement  or  condition. 
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9.0  List  of  Figures 


Figure  1  Schematic  illustrating  (i)  the  incident  wave,  Omo  originating  from 

the  source  at  position,  ps  (dashed  lines);  (ii)  the  first  order  scattered  wave, 
^scat  kl  arisinf>  from  the  first  heterogeneity  (solid  lines),  and  (iii)  the  second 

_ c\ 

order  wave,  k2  arising  from  re-scatter  of  the  first  order  wave  off  of  the 
second  heterogeneity. 

Figure  2  Schematic  detailing  the  geometry  used  in  the  calculation  of 
scattered  waves  from  analytical  expression.  The  centroid  of  the  cylinder  is 
the  origin  with  z  denoting  the  length,  angle  1)  denoting  the  angle  in  the  plane 
containing  the  source  and  detector,  and  r  denoting  the  radial  direction. 

Figure  3  Schematic  illustrating  the  geometry  of  the  experimental 
measurements. 

Figure  4  Experimental  values  of  0ki,k2  -  6inc  [degreesl  relative  to  the 
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igure  4  Experimental  values  of  0kl,k2  -  ©inc  [degrees]  relative  to  the 
absence  condition  as  a  function  of  distance  from  the  source-detector  pair. 
The  symbols  denote  individual  measurements  in  the  presence  of  two 
cylinders  separated  by  a  distance  (a)  6  mm,  (b)  10  mm,  and  (c)  20  mm 
while  the  line  connects  predictions  from  Eqn  (7)  and  measurements  of  0ki, 
0k2>  and  0inC.  The  error  bars  denote  the  propagation  of  measurement  errors 
(standard  deviation)  associated  with  0ki,  0k2,  and  ©inc-  The  x-axis  is 
reported  as  the  distance  between  the  wall  and  the  first  cylinder  (kl).  The 
astericks  denote  significant  difference  (p<0.005,  paired  student  t-test) 
between  the  values  experimentally  measured  and  obtained  from  Eqn  (7). 

Figure  5  Experimental  values  of  Mki,k2  /  Minc  [arbitrary  units]  relative  to 
the  absence  condition  as  a  function  of  distance  from  the  source-detector 
pair.  The  symbols  denote  individual  measurements  in  the  presence  of  two 
cylinders  separated  by  a  distance  (a)  6  mm,  (b)  10  mm,  and  (c)  20  mm 
while  the  line  connects  predictions  from  Eqn  (8)  and  measurements  of  Mki, 
Mk2,  and  Minc.  The  error  bars  denote  the  propagation  of  measurement 
errors  (standard  deviation)  associated  with  Mki,  Mk2,  and  Mine.  The  x-axis 
is  reported  as  the  distance  between  the  wall  and  the  first  cylinder  (kl). 
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Figure  6  Values  of  0ki,k2  -  ©inc  [degrees]  predicted  from  analytical  prediction 
of  first  order  (solid  line)  and  inclusive  of  second  order  (dashed  line)  scattering 
effects  as  a  function  of  distance  [cm]  from  the  source-detector  pair  for  two 
absorbing  cylinders  separated  by  (a)  6  mm,  (b)  10  mm,  and  (c)  20  mm.  The 
x-axis  is  reported  as  the  distance  between  the  wall  and  the  first  cylinder 
(kl)  and  the  phase-shift  reported  relative  to  the  absence  case. 

Figure  7  Finite  element  computations  of  0ki,k2  -  0inc  [degrees]  relative  to  an 
absence  condition  for  considering  only  first  order  (solid  line)  and  including 
second  order  (dashed  line)  perturbations  as  a  function  of  distance  [cm]  from 
the  source-detector  pair  for  absorbing  cylinders  separated  by  (a)  6  mm,  (b) 
10  mm  and  (c)  20  mm.  The  x-axis  is  reported  as  the  distance  between  the 
wall  and  the  first  cylinder  (kl)  and  the  phase-shift  reported  relative  to  the 
absence  case. 

Figure  8  Finite  element  computations  of  Mki,k2/Mjnc  [arbitrary  units] 
referenced  to  an  absence  condition  for  both  the  first  order  (solid  line)  and  the 
second  order  (dashed  line)  perturbations  as  a  function  of  distance  [cm]  from 
the  source-detector  pair  for  absorbing  cylinders  separated  by  (a)  6  mm,  (b) 
10  mm  and  (c)  20  mm.  The  x-axis  is  reported  as  the  distance  between  the 
wall  and  the  first  cylinder  (kl)  and  the  modulation  ratio  reported  relative  to 
the  absence  case. 
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